Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Phân tích số liệu và biểu đồ bằng Nguyễn Văn Tuấn Garvan Institute of Medical Research Sydney, Australia 1 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Mục lục 1 Tải R xuống và cài đặt vào máy tính 4 2 Tải R package và cài đặt vào máy tính 6 3 3.1 3.2 “Văn phạm” R Cách đặt tên trong R Hỗ trợ trong R 7 9 9 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Cách nhập dữ liệu vào R Nhập số liệu trực tiếp: c() Nhập số liệu trực tiếp: edit(data.frame()) Nhập số liệu từ một text file: read.table Nhập số liệu từ Excel Nhập số liệu từ SPSS Thông tin về số liệu Tạo dãy số bằng hàm seq, rep và gl 10 10 12 13 14 15 16 17 5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 Biên tập số liệu Tách rời số liệu: subset Chiết số liệu từ một data .frame Nhập hai data.frame thành một: merge Biến đổi số liệu (data coding) Biến đổi số liệu bằng cách dùng replace Biến đổi thành yếu tố (factor) Phân nhóm số liệu bằng cut2 (Hmisc) 19 19 20 21 22 23 23 24 6 6.1 6.2 Sử dụng R cho tính toán đơn giản Tính toán đơn giản Sử dụng R cho các phép tính ma trận 24 24 26 7 7.1 7.2 7.3 7.3.1 7.3.2 7.3.3 7.3.4 7.4 Sử dụng R cho tính toán xác suất Phép hoán vị (permutation) Biến số ngẫu nhiên và hàm phân phối Biến số ngẫu nhiên và hàm phân phối Hàm phân phối nhị phân (Binomial distribution) Hàm phân phối Poisson (Poisson distribution) Hàm phân phối chuẩn (Normal distribution) Hàm phân phối chuẩn chuẩn hóa (Standardized Normal distribution) Chọn mẫu ngẫu nhiên (random sampling) 31 31 32 32 33 35 36 38 41 8 8.1 8.2 8.3 8.4 8.5 8.5.1 8.5.2 8.6 8.7 8.7.1 8.8 Biểu đồ Số liệu cho phân tích biểu đồ Biểu đồ cho một biến số rời rạc (discrete variable): barplot Biểu đồ cho hai biến số rời rạc (discrete variable): barplot Biểu đồ hình tròn Biểu đồ cho một biến số liên tục: stripchart và hist Stripchart Histogram Biểu đồ hộp (boxplot) Phân tích biểu đồ cho hai biến liên tục Biểu đồ tán xạ (scatter plot) Phân tích Biểu đồ cho nhiều biến: pairs 42 42 44 45 46 47 47 48 49 50 50 53 2 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 8.9 Biểu đồ với sai số chuẩn (standard error) 54 9 9.1 9.2 9.3 9.3.1 9.3.2 9.4 9.5 9.6 9.7 9.8 9.9 9.10 9.10.1 9.10.2 Phân tích thống kê mô tả Thống kê mô tả (descriptive statistics, summary) Thống kê mô tả theo từng nhóm Kiểm định t (t.test) Kiểm định t một mẫu Kiểm định t hai mẫu Kiểm định Wilcoxon cho hai mẫu (wilcox.test) Kiểm định t cho các biến số theo cặp (paired t-test, t.test) Kiểm định Wilcoxon cho các biến số theo cặp (wilcox.test) Tần số (frequency) Kiểm định tỉ lệ (proportion test, prop.test, binom.test) So sánh hai tỉ lệ (prop.test, binom.test) So sánh nhiều tỉ lệ (prop.test, chisq.test) Kiểm định Chi bình ph ơng (Chi squared test, chisq.test) Kiểm định Fisher (Fisher’s exact test, fisher.test) 55 55 60 61 61 62 63 64 65 66 67 68 69 70 71 10 10.1 10.1.1 10.1.2 10.1.3 10.2 10.3 Phân tích hồi qui tuyến tính Hệ số t ơng quan Hệ số t ơng quan Pearson Hệ số t ơng quan Spearman Hệ số t ơng quan Kendall Mô hình của hồi qui tuyến tính đơn giản Mô hình hồi qui tuyến tính đa biến (multiple linear regression) 71 73 73 74 74 75 82 11 11.1 11.2 11.3 11.4 Phân tích ph ơng sai Phân tích ph ơng sai đơn giản (one-way analysis of variance) So sánh nhiều nhóm và điều chỉnh trị số p Phân tích bằng ph ơng pháp phi tham số Phân tích ph ơng sai hai chiều (two-way ANOVA) 85 85 87 90 91 12 12.1 12.2 12.3 Phân tích hồi qui logistic Mô hình hồi qui logistic Phân tích hồi qui logistic bằng R ớc tính xác suất bằng R 94 95 97 101 13 13.1 13.2 13.4 13.4.1 13.4.2 13.4.3 13.4.4 13.4.5 ớc tính cỡ mẫu (sample size estimation) Khái niệm về “power” Số liệu để ớc tính cỡ mẫu ớc tính cỡ mẫu ớc tính cỡ mẫu cho một chỉ số trung bình ớc tính cỡ mẫu cho so sánh hai số trung bình ớc tính cỡ mẫu cho phân tích ph ơng sai ớc tính cỡ mẫu để ớc tính một tỉ lệ ớc tính cỡ mẫu cho so sánh hai tỉ lệ 103 104 106 107 107 108 110 111 112 14 Tài liệu tham khảo 115 15 Thuật ngữ dùng trong sách 117 3 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Giới thiệu R Phân tích số liệu và biểu đồ th ng đ ợc tiến hành bằng các phần mềm thông dụng nh SAS, SPSS, Stata, Statistica, và S-Plus. Đây là những phần mềm đ ợc các công ti phần mềm phát triển và giới thiệu trên thị tr ng kho ng ba thập niên qua, và đã đ ợc các tr ng đ i học, các trung tâm nghiên cứu và công ti kĩ nghệ trên toàn thế giới sử dụng cho gi ng d y và nghiên cứu. Nh ng vì chi phí để sử dụng các phần mềm này tuơng đối đắt tiền (có khi lên đến hàng trăm ngàn đô-la mỗi năm), một số tr ng đ i học các n ớc đang phát triển (và ngay c một số n ớc đã phát triển) không có kh năng tài chính để sử dụng chúng một cách lâu dài. Do đó, các nhà nghiên cứu thống kê trên thế giới đã hợp tác với nhau để phát triển một phần mềm mới, với chủ tr ơng mã nguồn m , sao cho tất c các thành viên trong ngành thống kê học và toán học trên thế giới có thể sử dụng một cách thống nhất và hoàn toàn miễn phí. Năm 1996, trong một bài báo quan trọng về tính toán thống kê, hai nhà thống kê học Ross Ihaka và Robert Gentleman [lúc đó] thuộc Tr ng đ i học Auckland, New Zealand phát ho một ngôn ngữ mới cho phân tích thống kê mà họ đặt tên là R [1]. Sáng kiến này đ ợc rất nhiều nhà thống kê học trên thế giới tán thành và tham gia vào việc phát triển R. Cho đến nay, qua ch a đầy 10 năm phát triển, càng ngày càng có nhiều nhà thống kê học, toán học, nghiên cứu trong mọi lĩnh vực đã chuyển sang sử dụng R để phân tích dữ liệu khoa học. Trên toàn cầu, đã có một m ng l ới hơn một triệu ng i sử dụng R, và con số này đang tăng rất nhanh. Có thể nói trong vòng 10 năm nữa, vai trò của các phần mềm thống kê th ơng m i sẽ không còn lớn nh trong th i gian qua nữa. Vậy R là gì? Nói một cách ngắn gọn, R là một phần mềm sử dụng cho phân tích thống kê và vẽ biểu đồ. Thật ra, về b n chất, R là ngôn ngữ máy tính đa năng, có thể sử dụng cho nhiều mục tiêu khác nhau, từ tính toán đơn gi n, toán học gi i trí (recreational mathematics), tính toán ma trận (matrix), đến các phân tích thống kê phức t p. Vì là một ngôn ngữ, cho nên ng i ta có thể sử dụng R để phát triển thành các phần mềm chuyên môn cho một vấn đề tính toán cá biệt. Vì thế, những ai làm nghiên cứu khoa học, nhất là các n ớc còn nghèo khó nh n ớc ta, cần ph i học cách sử dụng R cho phân tích thống kê và đồ thị. Bài viết ngắn này sẽ h ớng dẫn b n đọc cách sử dụng R. Tôi gi định rằng b n đọc không biết gì về R, nh ng tôi kì vọng b n đọc biết qua về cách sử dụng máy tính. 1. Tải R xuống và cài đặt vào máy tính Để sử dụng R, việc đầu tiên là chúng ta ph i cài đặt R trong máy tính của mình. Để làm việc này, ta ph i truy nhập vào m ng và vào website có tên là “Comprehensive R Archive Network” (CRAN) sau đây: http://cran.R-project.org. 4 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Tài liệu cần t i về, tùy theo phiên b n, nh ng th ng có tên bắt đầu bằng mẫu tự R và số phiên b n (version). Chẳng h n nh phiên b n tôi sử dụng vào cuối năm 2005 là 2.2.1, nên tên của tài liệu cần t i là: R-2.2.1-win32.zip Tài liệu này kho ng 26 MB, và địa chỉ cụ thể để t i là: http://cran.r-project.org/bin/windows/base/R-2.2.1-win32.exe T i website này, chúng ta có thể tìm thấy rất nhiều tài liệu chỉ dẫn cách sử dụng R, đủ trình độ, từ sơ đẳng đến cao cấp. Nếu ch a quen với tiếng Anh, tài liệu này của tôi có thể cung cấp những thông tin cần thiết để sử dụng mà không cần ph i đọc các tài liệu khác. Khi đã t i R xuống máy tính, b ớc kế tiếp là cài đặt (set-up) vào máy tính. Để làm việc này, chúng ta chỉ đơn gi n nhấn chuột vào tài liệu trên và làm theo h ớng dẫn cách cài đặt trên màn hình. Đây là một b ớc rất đơn gi n, chỉ cần 1 phút là việc cài đặt R có thể hoàn tất. Sau khi hoàn tất việc cài đặt, một icon R 2.2.1.lnk sẽ xuất hiện trên desktop của máy tính. Đến đây thì chúng ta đã sẵn sàng sử dụng R. Có thể nhấp chuột vào icon này và chúng ta sẽ có một window nh sau: 5 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 2. Tải R package và cài đặt vào máy tính R cung cấp cho chúng ta một “ngôn ngữ” máy tính và một số function để làm các phân tích căn b n và đơn gi n. Nếu muốn làm những phân tích phức t p hơn, chúng ta cần ph i t i về máy tính một số package khác. Package là một phần mềm nhỏ đ ợc các nhà thống kê phát triển để gi i quyết một vấn đề cụ thể, và có thể ch y trong hệ thống R. Chẳng h n nh để phân tích hồi qui tuyến tính, R có function lm để sử dụng cho mục đích này, nh ng để làm các phân tích sâu hơn và phức t p hơn, chúng ta cần đến các package nh lme4. Các package này cần ph i đ ợc t i về và cài đặt vào máy tính. Địa chỉ để t i các package vẫn là: http://cran.r-project.org, rồi bấm vào phần “Packages” xuất hiện bên trái của mục lục trang web. Theo tôi, một số package cần t i về máy tính để sử dụng cho các phân tích dịch tễ học là: Tên package trellis lattice Hmisc Design Epi epitools Foreign Rmeta meta Ch c năng Dùng để vẽ đồ thị và làm cho đồ thị đẹp hơn Dùng để vẽ đồ thị và làm cho đồ thị đẹp hơn Một số ph ơng pháp mô hình dữ liệu của F. Harrell Một số mô hình thiết kế nghiên cứu của F. Harrell Dùng cho các phân tích dịch tễ học Một package khác chuyên cho các phân tích dịch tễ học Dùng để nhập dữ liệu từ các phần mềm khác nh SPSS, Stata, SAS, v.v… Dùng cho phân tích tổng hợp (meta-analysis) Một package khác cho phân tích tổng hợp 6 Phân tích số liệu và biểu đồ bằng R survival Nguyễn Văn Tuấn Chuyên dùng cho phân tích theo mô hình Cox (Cox’s proportional hazard model) Package dùng cho các phân tích thống kê trong lĩnh vực xã hội học Package dùng cho phân tích số liệu di truyền học Bayesian Model Average Zelig Genetics BMA Các package này có thể cài đặt trực tuyến bằng cách chọn Install packages trong phần packages của R nh hình d ới đây. Ngoài ra, nếu package đã đ ợc t i xuống máy tính cá nhân, việc cài đặt có thể nhanh hơn bằng cách chọn Install package(s) from local zip file cũng trong phần packages (xem hình d ới đây). 3. “Văn phạm” R R là một ngôn ngữ t ơng tác (interactive language), có nghĩa là khi chúng ta ra lệnh, và nếu lệnh theo đúng “văn ph m”, R sẽ “đáp” l i bằng một kết qu . Và, sự t ơng tác tiếp tục cho đến khi chúng ta đ t đ ợc yêu cầu. “Văn ph m” chung của R là một lệnh (command) hay function (tôi sẽ thỉnh tho ng đề cập đến là “hàm”). Mà đã là hàm thì ph i có thông số; cho nên theo sau hàm là những thông số mà chúng ta ph i cung cấp. Cú pháp chung của R là nh sau: đối tượng <- hàm(thông số 1, thông số 2, …, thông số n) 7 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Chẳng h n nh : > reg <- lm(y ~ x) thì reg là một đối t ợng (object), còn lm là một hàm, và y ~ x là thông số của hàm. Hay: > setwd(“c:/works/stats”) thì setwd là một hàm, còn “c:/works/stats” là thông số của hàm. Để biết một hàm cần có những thông số nào, chúng ta dùng lệnh args(x), (args viết tắt chữ arguments) mà trong đó x là một hàm chúng ta cần biết: > args(lm) function (formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...) NULL R là một ngôn ngữ “đối t ợng” (object oriented language). Điều này có nghĩa là các dữ liệu trong R đ ợc chứa trong object. Định h ớng này cũng có vài nh h ng đến cách viết của R. Chẳng h n nh thay vì viết x = 5 nh thông th ng chúng ta vẫn viết, thì R yêu cầu viết là x == 5. Đối với R, x = 5 t ơng đ ơng với x <- 5. Cách viết sau (dùng kí hiệu <-) đ ợc khuyến khích hơn là cách viết tr ớc (=). Chẳng h n nh : > x <- rnorm(10) có nghĩa là mô phỏng 10 số liệu và chứa trong object x. Chúng ta cũng có thể viết x = rnorm(10). Một số kí hiệu hay dùng trong R là: x == 5 x != 5 y < x x > y z <= 7 p >= 1 is.na(x) A & B A | B ! x bằng 5 x không bằng 5 y nhỏ hơn x x lớn hơn y z nhỏ hơn hoặc bằng 7 p lớn hơn hoặc bằng 1 Có ph i x là biến số trống không (missing value) A và B (AND) A hoặc B (OR) Không là (NOT) 8 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Với R, tất c các câu chữ hay lệnh sau kí hiệu # đều không có hiệu ứng, vì # là kí hiệu dành cho ng i sử dụng thêm vào các ghi chú, ví dụ: > # lệnh sau đây sẽ mô phỏng 10 giá trị normal > x <- rnorm(10) 3.1 Cách đặt tên trong R Đặt tên một đối t ợng (object) hay một biến số (variable) trong R khá linh ho t, vì R không có nhiều giới h n nh các phần mềm khác. Tên một object ph i đ ợc viết liền nhau (tức không đ ợc cách r i bằng một kho ng trống). Chẳng h n nh R chấp nhận myobject nh ng không chấp nhận my object. > myobject <- rnorm(10) > my object <- rnorm(10) Error: syntax error in "my object" Nh ng đôi khi tên myobject khó đọc, cho nên chúng ta nên tác r i bằng “.” Nh my.object. > my.object <- rnorm(10) Một điều quan trọng cần l u ý là R phân biệt mẫu tự viết hoa và viết th My.object khác với my.object. Ví dụ: ng. Cho nên > My.object.u <- 15 > my.object.L <- 5 > My.object.u + my.object.L [1] 20 Một vài điều cần l u ý khi đặt tên trong R là: • • Không nên đặt tên một biến số hay variable bằng kí hiệu “_” (underscore) nh my_object hay my-object. Không nên đặt tên một object giống nh một biến số trong một dữ liệu. Ví dụ, nếu chúng ta có một data.frame (dữ liệu hay dataset) với biến số age trong đó, thì không nên có một object trùng tên age, tức là không nên viết: age <age. Tuy nhiên, nếu data.frame tên là data thì chúng ta có thể đề cập đến biến số age với một kí tự $ nh sau: data$age. (Tức là biến số age trong data.frame data), và trong tr ng hợp đó, age <- data$age có thể chấp nhận đ ợc. 3.2 Hỗ trợ trong R 9 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ngoài lệnh args() R còn cung cấp lệnh help() để ng i sử dụng có thể hiểu “văn ph m” của từng hàm. Chẳng h n nh muốn biết hàm lm có những thông số (arguments) nào, chúng ta chỉ đơn gi n lệnh: > help(lm) hay > ?lm Một cửa sổ sẽ hiện ra bên ph i của màn hình chỉ rõ cách sử dụng ra sao và thậm chí có c ví dụ. B n đọc có thể đơn gi n copy và dán ví dụ vào R để xem cách vận hành. Tr ớc khi sử dụng R, ngoài sách này nếu cần b n đọc có thể đọc qua phần chỉ dẫn có sẵn trong R bằng cách chọn mục help và sau đó chọn Html help nh hình d ới đây để biết thêm chi tiết. B n đọc cũng có thể copy và dán các lệnh trong mục này vào R để xem cho biết cách vận hành của R. 4. Cách nhập dữ liệu vào R Muốn làm phân tích dữ liệu bằng R, chúng ta ph i có sẵn dữ liệu d ng mà R có thể hiểu đ ợc để xử lí. Dữ liệu mà R hiểu đ ợc ph i là dữ liệu trong một data.frame. Có nhiều cách để nhập số liệu vào một data.frame trong R, từ nhập trực tiếp đến nhập từ các nguồn khác nhau. Sau đây là những cách thông dụng nhất: 4.1 Nhập số liệu trực tiếp: c() Ví dụ 1: chúng ta có số liệu về độ tuổi và insulin cho 10 bệnh nhân nh sau, và muốn nhập vào R. 50 62 60 40 48 47 57 70 48 67 16.5 10.8 32.3 19.3 14.2 11.3 15.5 15.8 16.2 11.2 Chúng ta có thể sử dụng function có tên c nh sau: > age <- c(50,62, 60,40,48,47,57,70,48,67) > insulin <- c(16.5,10.8,32.3,19.3,14.2,11.3,15.5,15.8,16.2,11.2) 10 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Lệnh thứ nhất cho R biết rằng chúng ta muốn t o ra một cột dữ liệu (từ nay tôi sẽ gọi là biến số, tức variable) có tên là age, và lệnh thứ hai là t o ra một cột khác có tên là insulin. Tất nhiên, chúng ta có thể lấy một tên khác mà mình thích. Chúng ta dùng function c (viết tắt của chữ concatenation – có nghĩa là “móc nối vào nhau”) để nhập dữ liệu. Chú ý rằng mỗi số liệu cho mỗi bệnh nhân đ ợc cách nhau bằng một dấu phẩy. Kí hiệu insulin <- (cũng có thể viết là insulin =) có nghĩa là các số liệu theo sau sẽ có nằm trong biến số insulin. Chúng ta sẽ gặp kí hiệu này rất nhiều lần trong khi sử dụng R. R là một ngôn ngữ cấu trúc theo d ng đối t ợng (thuật ngữ chuyên môn là “object-oriented language”), vì mỗi cột số liệu hay mỗi một data.frame là một đối t ợng (object) đối với R. Vì thế, age và insulin là hai đối t ợng riêng lẻ. Bây gi chúng ta cần ph i nhập hai đối t ợng này thành một data.frame để R có thể xử lí sau này. Để làm việc này chúng ta cần đến function data.frame: > tuan <- data.frame(age, insulin) Trong lệnh này, chúng ta muốn cho R biết rằng nhập hai cột (hay hai đối t ợng) age và insulin vào một đối t ợng có tên là tuan. Đến đây thì chúng ta đã có một đối t ợng hoàn chỉnh để tiến hành phân tích thống kê. Để kiểm tra xem trong tuan có gì, chúng ta chỉ cần đơn gi n gõ: > tuan Và R sẽ báo cáo: 1 2 3 4 5 6 7 8 9 10 age insulin 50 16.5 62 10.8 60 32.3 40 19.3 48 14.2 47 11.3 57 15.5 70 15.8 48 16.2 67 11.2 Nếu chúng ta muốn l u l i các số liệu này trong một file theo d ng R, chúng ta cần dùng lệnh save. Gi dụ nh chúng ta muốn l u số liệu trong directory có tên là “c:\works\insulin”, chúng ta cần gõ nh sau: > setwd(“c:/works/insulin”) > save(tuan, file=”tuan.rda”) 11 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Lệnh đầu tiên (setwd – chữ wd có nghĩa là working directory) cho R biết rằng chúng ta muốn l u các số liệu trong directory có tên là “c:\works\insulin”. L u ý rằng thông th ng Windows dùng dấu backward slash “/”, nh ng trong R chúng ta dùng dấu forward slash “/”. Lệnh thứ hai (save) cho R biết rằng các số liệu trong đối t ợng tuan sẽ l u trong file có tên là “tuan.rda”). Sau khi gõ xong hai lệnh trên, một file có tên tuan.rda sẽ có mặt trong directory đó. 4.2 Nhập số liệu trực tiếp: edit(data.frame()) Ví dụ 1 (tiếp tục): chúng ta có thể nhập số liệu về độ tuổi và insulin cho 10 bệnh nhân bằng một function rất có ích, đó là: edit(data.frame()). Với function này, R sẽ cung cấp cho chúng ta một window mới với một dãy cột và dòng giống nh Excel, và chúng ta có thể nhập số liệu trong b ng đó. Ví dụ: > ins <- edit(data.frame()) Chúng ta sẽ có một cửa sổ nh sau: đây, R không biết chúng ta có biến số nào, cho nên R liệt kê các biến số var1, var2, v.v… Nhấp chuột vào cột var1 và thay đổi bằng cách gõ vào đó age. Nhấp chuột vào cột var2 và thay đổi bằng cách gõ vào đó insulin. Sau đó gõ số liệu cho 12 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn từng cột. Sau khi xong, bấm nút chéo X góc ph i của spreadsheet, chúng ta sẽ có một data.frame tên ins với hai biến số age và insulin. 4.3 Nhập số liệu từ một text file: read.table Ví dụ 2: Chúng ta thu thập số liệu về độ tuổi và cholesterol từ một nghiên cứu 50 bệnh nhân mắc bệnh cao huyết áp. Các số liệu này đ ợc l u trong một text file có tên là chol.txt t i directory c:\works\insulin. Số liệu này nh sau: cột 1 là mã số của bệnh nhân, cột 2 là giới tính, cột 3 là body mass index (bmi), cột 4 là HDL cholesterol (viết tắt là hdl), kế đến là LDL cholesterol, total cholesterol (tc) và triglycerides (tg). id 1 2 3 4 5 6 7 8 9 10 ... 46 47 48 49 50 sex Nam Nu Nu Nam Nam Nu Nam Nam Nam Nu age 57 64 60 65 47 65 76 61 59 57 bmi 17 18 18 18 18 18 19 19 19 19 hdl 5.000 4.380 3.360 5.920 6.250 4.150 0.737 7.170 6.942 5.000 ldl 2.0 3.0 3.0 4.0 2.1 3.0 3.0 3.0 3.0 2.0 tc 4.0 3.5 4.7 7.7 5.0 4.2 5.9 6.1 5.9 4.0 tg 1.1 2.1 0.8 1.1 2.1 1.5 2.6 1.5 5.4 1.9 Nu Nam Nam Nu Nu 52 64 45 64 62 24 24 24 25 25 3.360 7.170 7.880 7.360 7.750 2.0 1.0 4.0 4.6 4.0 3.7 6.1 6.7 8.1 6.2 1.2 1.9 3.3 4.0 2.5 Chúng ta muốn nhập các dữ liệu này vào R để tiện việc phân tích sau này. Chúng ta sẽ sử dụng lệnh read.table nh sau: > setwd(“c:/works/insulin”) > chol <- read.table("chol.txt", header=TRUE) Lệnh thứ nhất chúng ta muốn đ m b o R truy nhập đúng directory mà số liệu đang đ ợc l u giữ. Lệnh thứ hai yêu cầu R nhập số liệu từ file có tên là “chol.txt” (trong directory c:\works\insulin) và cho vào đối t ợng chol. Trong lệnh này, header=TRUE có nghĩa là yêu cầu R đọc dòng đầu tiên trong file đó nh là tên của từng cột dữ kiện. Chúng ta có thể kiểm tra xem R đã đọc hết các dữ liệu hay ch a bằng cách ra lệnh: > chol Hay 13 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > names(chol) R sẽ cho biết có các cột nh sau trong dữ liệu (names là lệnh hỏi trong dữ liệu có những cột nào và tên gì): [1] "id" "sex" "age" "bmi" "hdl" "ldl" "tc" "tg" Bây gi chúng ta có thể l u dữ liệu d ới d ng R để xử lí sau này bằng cách ra lệnh: > save(chol, file="chol.rda") 4.4 Nhập số liệu từ Excel: read.csv Để nhập số liệu từ phần mềm Excel, chúng ta cần tiến hành 2 b ớc: • • B ớc 1: Dùng lệnh “Save as” trong Excel và l u số liệu d ới d ng “csv”; B ớc 2: Dùng R (lệnh read.csv) để nhập dữ liệu d ng csv. Ví dụ 3: Một dữ liệu gồm các cột sau đây đang đ ợc l u trong Excel, và chúng ta muốn chuyển vào R để phân tích. Dữ liệu này có tên là excel.xls. ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Age 18 28 20 21 28 23 20 20 20 20 22 27 26 33 34 32 28 18 26 27 Sex 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 Ethnicity 1 1 1 1 1 4 1 1 1 1 1 2 1 1 3 1 1 2 2 2 IGFI 148.27 114.50 109.82 112.13 102.86 129.59 142.50 118.69 197.69 163.69 144.81 141.60 161.80 89.20 161.80 148.50 157.70 222.90 186.70 167.56 IGFBP3 5.14 5.23 4.33 4.38 4.04 4.16 3.85 3.44 4.12 3.96 3.63 3.48 4.10 2.82 3.80 3.72 3.98 3.98 4.64 3.56 ALS 316.00 296.42 269.82 247.96 240.04 266.95 300.86 277.46 335.23 306.83 295.46 231.20 244.80 177.20 243.60 234.80 224.80 281.40 340.80 321.12 PINP 61.84 98.64 93.26 101.59 58.77 48.93 135.62 79.51 57.25 74.03 68.26 56.78 75.75 48.57 50.68 83.98 60.42 74.17 38.05 30.18 ICTP 5.81 4.96 7.74 6.66 4.62 5.32 8.78 7.19 6.21 4.95 4.54 4.47 6.27 3.58 3.52 4.85 4.89 6.43 5.12 4.78 P3NP 4.21 5.33 4.56 4.61 4.95 3.82 6.75 5.11 4.44 4.84 3.70 4.07 5.26 3.68 3.35 3.80 4.09 5.84 5.77 6.12 Việc đầu tiên là chúng ta cần làm, nh nói trên, là vào Excel để l u d ới d ng csv: • Vào Excel, chọn File Æ Save as • Chọn Save as type “CSV (Comma delimited)” 14 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Sau khi xong, chúng ta sẽ có một file với tên “excel.csv” trong directory “c:\works\insulin”. Việc thứ hai là vào R và ra những lệnh sau đây: > setwd(“c:/works/insulin”) > gh <- read.csv ("excel.txt", header=TRUE) Lệnh thứ hai read.csv yêu cầu R đọc số liệu từ “excel.csv”, dùng dòng thứ nhất là tên cột, và l u các số liệu này trong một object có tên là gh. Bây gi chúng ta có thể l u gh d ới d ng R để xử lí sau này bằng lệnh sau đây: > save(gh, file="gh.rda") 4.5 Nhập số liệu từ một SPSS: read.spss Phần mềm thống kê SPSS l u dữ liệu d ới d ng “sav”. Chẳng h n nh nếu chúng ta đã có một dữ liệu có tên là testo.sav trong directory c:\works\insulin, và muốn chuyển dữ liệu này sang d ng R có thể hiểu đ ợc, chúng ta cần sử dụng lệnh read.spss trong package có tên là foreign. Các lệnh sau đây sẽ hoàn tất dễ dàng việc này: Việc đầu tiên chúng ta cho truy nhập foreign bằng lệnh library: 15 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > library(foreign) Việc thứ hai là lệnh read.spss: > setwd(“c:/works/insulin”) > testo <- read.spss(“testo.sav”, to.data.frame=TRUE) Lệnh thứ hai read.spss yêu cầu R đọc số liệu từ “testo.sav”, và cho vào một data.frame có tên là testo. Bây gi chúng ta có thể l u testo d ới d ng R để xử lí sau này bằng lệnh sau đây: > save(testo, file="testo.rda") 4.6 Thông tin về dữ liệu Gi dụ nh chúng ta đã nhập số liệu vào một data.frame có tên là chol nh trong ví dụ 1. Để tìm hiểu xem trong dữ liệu này có gì, chúng ta có thể nhập vào R nh sau: • Dẫn cho R biết chúng ta muốn xử lí chol bằng cách dùng lệnh attach(arg) với arg là tên của dữ liệu.. > attach(chol) • Chúng ta có thể kiểm tra xem chol có ph i là một data.frame không bằng lệnh is.data.frame(arg) với arg là tên của dữ liệu. Ví dụ: > is.data.frame(chol) [1] TRUE R cho biết chol qu là một data.frame. • Có bao nhiêu cột (hay variable = biến số) và dòng số liệu (observations) trong dữ liệu này? Chúng ta dùng lệnh dim(arg) với arg là tên của dữ liệu. (dim viết tắt chữ dimension). Ví dụ (kết qu của R trình bày ngay sau khi chúng ta gõ lệnh): > dim(chol) [1] 50 8 • Nh vậy, chúng ta có 50 dòng và 8 cột (hay biến số). Vậy những biến số này tên gì? Chúng ta dùng lệnh names(arg) với arg là tên của dữ liệu. Ví dụ: > names(chol) [1] "id" "sex" "age" "bmi" "hdl" "ldl" "tc" 16 "tg" Phân tích số liệu và biểu đồ bằng R • Nguyễn Văn Tuấn Trong biến số sex, chúng ta có bao nhiêu nam và nữ? Để tr l i câu hỏi này, chúng ta có thể dùng lệnh table(arg) với arg là tên của biến số. Ví dụ: > table(sex) sex nam Nam 1 21 Nu 28 Kết qu cho thấy dữ liệu này có 21 nam và 28 nữ. 4.7 Tạo dãy số bằng hàm seq, rep và gl R còn có công dụng t o ra những dãy số rất tiện cho việc mô phỏng và thiết kế thí nghiệm. Những hàm thông th ng cho dãy số là seq (sequence), rep (repetition) và gl (generating levels): Áp dụng seq • T o ra một vector số từ 1 đến 12: > x <- (1:12) > x [1] 1 2 3 > seq(12) [1] 1 2 • 3 4 4 5 5 6 6 7 7 8 8 9 10 11 12 9 10 11 12 T o ra một vector số từ 12 đến 5: > x <- (12:5) > x [1] 12 11 10 9 8 7 > seq(12,7) [1] 12 11 10 8 7 9 6 5 Công thức chung của hàm seq là seq(from, to, by= ) hay seq(from, to, length.out= ). Cách sử dụng sẽ đ ợc minh ho bằng vài ví dụ sau đây: • T o ra một vector số từ 4 đến 6 với kho ng cách bằng 0.25: • T o ra một vector 10 số, với số nhỏ nhất là 2 và số lớn nhất là 15 > seq(4, 6, 0.25) [1] 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 > seq(length=10, from=2, to=15) [1] 2.000000 3.444444 4.888889 6.333333 10.666667 12.111111 13.555556 15.000000 17 7.777778 9.222222 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Áp dụng rep Công thức của hàm rep là rep(x, times, ...), trong đó, x là một biến số và times là số lần lặp l i. Ví dụ: • T o ra số 10, 3 lần: • T o ra số 1 đến 4, 3 lần: • T o ra số 1.2, 2.7, 4.8, 5 lần: • T o ra số 1.2, 2.7, 4.8, 5 lần: > rep(10, 3) [1] 10 10 10 > rep(c(1:4), 3) [1] 1 2 3 4 1 2 3 4 1 2 3 4 > rep(c(1.2, 2.7, 4.8), 5) [1] 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 > rep(c(1.2, 2.7, 4.8), 5) [1] 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 1.2 2.7 4.8 Áp dụng gl gl đ ợc áp dụng để t o ra một biến thứ bậc (categorical variable), tức biến không để tính toán, mà là đếm. Công thức chung của hàm gl là gl(n, k, length = n*k, labels = 1:n, ordered = FALSE) và cách sử dụng sẽ đ ợc minh ho bằng vài ví dụ sau đây: • T o ra biến gồm bậc 1 và 2; mỗi bậc đ ợc lặp l i 8 lần: > gl(2, 8) [1] 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 Levels: 1 2 Hay một biến gồm bậc 1, 2 và 3; mỗi bậc đ ợc lặp l i 5 lần: > gl(3, 5) [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 Levels: 1 2 3 • T o ra biến gồm bậc 1 và 2; mỗi bậc đ ợc lặp l i 10 lần (do đó length=20): > gl(2, 10, length=20) [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 Levels: 1 2 Hay: > gl(2, 2, length=20) [1] 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 Levels: 1 2 • Cho thêm kí hiệu: 18 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > gl(2, 5, label=c("C", "T")) [1] C C C C C T T T T T Levels: C T • T o một biến gồm 4 bậc 1, 2, 3, 4. Mỗi bậc lặp l i 2 lần. > rep(1:4, c(2,2,2,2)) [1] 1 1 2 2 3 3 4 4 Cũng t ơng đ ơng với: > rep(1:4, each = 2) [1] 1 1 2 2 3 3 4 4 • Với ngày gi tháng: > x <- .leap.seconds[1:3] > rep(x, 2) [1] "1972-06-30 17:00:00 Pacific Standard Time" "1972-12-31 16:00:00 Pacific Standard Time" [3] "1973-12-31 16:00:00 Pacific Standard Time" "1972-06-30 17:00:00 Pacific Standard Time" [5] "1972-12-31 16:00:00 Pacific Standard Time" "1973-12-31 16:00:00 Pacific Standard Time" > rep(as.POSIXlt(x), rep(2, 3)) [1] "1972-06-30 17:00:00 Pacific Standard Time" "1972-06-30 17:00:00 Pacific Standard Time" [3] "1972-12-31 16:00:00 Pacific Standard Time" "1972-12-31 16:00:00 Pacific Standard Time" [5] "1973-12-31 16:00:00 Pacific Standard Time" "1973-12-31 16:00:00 Pacific Standard Time" 5. Biên tập số liệu 5.1 Tách rời dữ liệu: subset Chúng ta sẽ quay l i với dữ liệu chol trong ví dụ 1. Để tiện việc theo dõi và hiểu “câu chuyện”, tôi xin nhắc l i rằng chứng ta đã nhập số liệu vào trong một dữ liệu R có tên là chol từ một text file có tên là chol.txt: > setwd(“c:/works/insulin”) > chol <- read.table(“chol.txt”, header=TRUE) > attach(chol) Nếu chúng ta, vì một lí do nào đó, chỉ muốn phân tích riêng cho nam giới, chúng ta có thể tách chol ra thành hai data.frame, t m gọi là nam và nu. Để làm chuyện này, chúng ta dùng lệnh subset(data, cond), trong đó data là data.frame mà chúng ta muốn tách r i, và cond là điều kiện. Ví dụ: > nam <- subset(chol, sex==”Nam”) > nu <- subset(chol, sex==”Nu”) 19 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Sau khi ra hai lệnh này, chúng ta đã có 2 dữ liệu (hai data.frame) mới tên là nam và nu. Chú ý điều kiện sex == “Nam” và sex == “Nu” chúng ta dùng == thay vì = để chỉ điều kiện chính xác. Tất nhiên, chúng ta cũng có thể tách dữ liệu thành nhiều data.frame khác nhau với những điều kiện dựa vào các biến số khác. Chẳng h n nh lệnh sau đây t o ra một data.frame mới tên là old với những bệnh nhân trên 60 tuổi: > old <- subset(chol, age>=60) > dim(old) [1] 25 8 Hay một data.frame mới với những bệnh nhân trên 60 tuổi và nam giới: > n60 <- subset(chol, age>=60 & sex==”Nam”) > dim(n60) [1] 9 8 5.2 Chiết số liệu từ một data .frame Trong chol có 8 biến số. Chúng ta có thể chiết dữ liệu chol và chỉ giữ l i những biến số cần thiết nh mã số (id), độ tuổi (age) và total cholestrol (tc). Để ý từ lệnh names(chol) rằng biến số id là cột số 1, age là cột số 3, và biến số tc là cột số 7. Chúng ta có thể dùng lệnh sau đây: > data2 <- chol[, c(1,3,7)] đây, chúng ta lệnh cho R biết rằng chúng ta muốn chọn cột số 1, 3 và 7, và đ a tất c số liệu của hai cột này vào data.frame mới có tên là data2. Chú ý chúng ta sử dụng ngoặc kép vuông [] chứ không ph i ngoặc kép vòng (), vì chol không ph i làm một function. Dấu phẩy phía tr ớc c, có nghĩa là chúng ta chọn tất c các dòng số liệu trong data.frame chol. Nh ng nếu chúng ta chỉ muốn chọn 10 dòng số liệu đầu tiên, thì lệnh sẽ là: > data3 <- chol[1:10, c(1,3,7)] > print(data3) 1 2 3 4 5 6 7 8 id 1 2 3 4 5 6 7 8 sex Nam Nu Nu Nam Nam Nu Nam Nam tc 4.0 3.5 4.7 7.7 5.0 4.2 5.9 6.1 20 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 9 9 Nam 5.9 10 10 Nu 4.0 Chú ý lệnh print(arg) đơn gi n liệt kê tất c số liệu trong data.frame arg. Thật ra, chúng ta chỉ cần đơn gi n gõ data3, kết qu cũng giống y nh print(data3). 5.3 Nhập hai data.frame thành một: merge Gi dụ nh chúng ta có dữ liệu chứa trong hai data.frame. Dữ liệu thứ nhất tên là d1 gồm 3 cột: id, sex, tc nh sau: id sex tc 1 Nam 4.0 2 Nu 3.5 3 Nu 4.7 4 Nam 7.7 5 Nam 5.0 6 Nu 4.2 7 Nam 5.9 8 Nam 6.1 9 Nam 5.9 10 Nu 4.0 Dữ liệu thứ hai tên là d2 gồm 3 cột: id, sex, tg nh sau: id 1 2 3 4 5 6 7 8 9 10 11 sex Nam Nu Nu Nam Nam Nu Nam Nam Nam Nu Nu tg 1.1 2.1 0.8 1.1 2.1 1.5 2.6 1.5 5.4 1.9 1.7 Hai dữ liệu này có chung hai biến số id và sex. Nh ng dữ liệu d1 có 10 dòng, còn dữ liệu d2 có 11 dòng. Chúng ta có thể nhập hai dữ liệu thành một data.frame bằng cách dùng lệnh merge nh sau: > d <- merge(d1, d2, by="id", all=TRUE) > d id sex.x tc sex.y tg 21 Phân tích số liệu và biểu đồ bằng R 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 Nam Nu Nu Nam Nam Nu Nam Nam Nam Nu <NA> 4.0 3.5 4.7 7.7 5.0 4.2 5.9 6.1 5.9 4.0 NA Nam Nu Nu Nam Nam Nu Nam Nam Nam Nu Nu Nguyễn Văn Tuấn 1.1 2.1 0.8 1.1 2.1 1.5 2.6 1.5 5.4 1.9 1.7 Trong lệnh merge, chúng ta yêu cầu R nhập 2 dữ liệu d1 và d2 thành một và đ a vào data.frame mới tên là d, và dùng biến số id làm chuẩn. Chúng ta để ý thấy bệnh nhân số 11 không có số liệu cho tc, cho nên R cho là NA (một d ng “not available”). 5.4 Biến đổi số liệu (data coding) Trong việc xử lí số liệu dịch tễ học, nhiều khi chúng ta cần ph i biến đổi số liệu từ biến liên tục sang biến mang tính cách phân lo i. Chẳng h n nh trong chẩn đoán loãng x ơng, những phụ nữ có chỉ số T của mật độ chất khoáng trong x ơng (bone mineral density hay BMD) bằng hay thấp hơn -2.5 đ ợc xem là “loãng x ơng”, những ai có BMD giữa -2.5 và -1.0 là “xốp x ơng” (osteopenia), và trên -1.0 là “bình th ng”. Ví dụ, chúng ta có số liệu BMD từ 10 bệnh nhân nh sau: -0.92, 0.21, 0.17, -3.21, -1.80, -2.60, -2.00, 1.71, 2.12, -2.11 Để nhập các số liệu này vào R chúng ta có thể sử dụng function c nh sau: bmd <- c(-0.92,0.21,0.17,-3.21,-1.80,-2.60,-2.00,1.71,2.12,-2.11) Để phân lo i 3 nhóm loãng x ơng, xốp x ơng, và bình th ng, chúng ta có thể dùng mã số 1, 2 và 3. Nói cách khác, chúng ta muốn t o nên một biến số khác (hãy gọi là diagnosis) gồm 3 giá trị trên dựa vào giá trị của bmd. Để làm việc này, chúng ta sử dụng lệnh: # tạm thời cho biến số diagnosis bằng bmd > diagnosis <- bmd # > > > biến đổi bmd thành diagnosis diagnosis[bmd <= -2.5] <- 1 diagnosis[bmd > -2.5 & bmd <= 1.0] <- 2 diagnosis[bmd > -1.0] <- 3 # tạo thành một data frame > data <- data.frame(bmd, diagnosis) # liệt kê để kiểm tra xem lệnh có hiệu quả không > data 22 Phân tích số liệu và biểu đồ bằng R 1 2 3 4 5 6 7 8 9 10 Nguyễn Văn Tuấn bmd diagnosis -0.92 3 0.21 3 0.17 3 -3.21 1 -1.80 2 -2.60 1 -2.00 2 1.71 3 2.12 3 -2.11 2 5.5 Biến đổi số liệu bằng cách dùng replace Một cách biến đổi số liệu khác là dùng replace, dù cách này có vẻ r m rà chút ít. Tiếp tục ví dụ trên, chúng ta biến đổi từ bmd sang diagnosis nh sau: > > > > diagnosis diagnosis diagnosis diagnosis <<<<- bmd replace(diagnosis, bmd <= -2.5, 1) replace(diagnosis, bmd > -2.5 & bmd <= 1.0, 2) replace(diagnosis, bmd > -1.0, 3) 5.6 Biến đổi thành yếu tố (factor) Trong phân tích thống kê, chúng ta phân biệt một biến số mang tính yếu tố (factor) và biến số liên tục bình th ng. Biến số yếu tố không thể dùng để tính toán nh cộng trừ nhân chia, nh ng biến số số học có thể sử dụng để tính toán. Chẳng h n nh trong ví dụ bmd và diagnosis trên, diagnosis là yếu tố vì giá trị trung bình giữa 1 và 2 chẳng có ý nghĩa thực tế gì c ; còn bmd là biến số số học. Nh ng hiện nay, diagnosis đ ợc xem là một biến số số học. Để biến thành biến số yếu tố, chúng ta cần sử dụng function factor nh sau: > diag <- factor(diagnosis) > diag [1] 3 3 3 1 2 1 2 3 3 2 Levels: 1 2 3 Chú ý R bây gi thông báo cho chúng ta biết diag có 3 bậc: 1, 2 và 3. Nếu chúng ta yêu cầu R tính số trung bình của diag, R sẽ không làm theo yêu cầu này, vì đó không ph i là một biến số số học: > mean(diag) [1] NA Warning message: argument is not numeric or logical: returning NA in: mean.default(diag) Dĩ nhiên, chúng ta có thể tính giá trị trung bình của diagnosis: 23 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > mean(diagnosis) [1] 2.3 nh ng kết qu 2.3 này không có ý nghĩa gì trong thực tế c . 5.7 Phân nhóm số liệu bằng cut2 (Hmisc) Trong phân tích thống kê, có khi chúng ta cần ph i phân chia một biến số liên tục thành nhiều nhóm dựa vào phân phối của biến số. Chẳng h n nh đối với biến số bmd chúng ta có thể “cắt” dãy số thành 3 nhóm t ơng đ ơng nhau bằng cách dùng function cut2 (trong th viện Hmisc) nh sau: > # nhập thư viện Hmisc để có thể dùng function cut2 > library(Hmisc) > bmd <- c(-0.92,0.21,0.17,-3.21,-1.80,-2.60,-2.00,1.71,2.12,-2.11) > # chia biến số bmd thành 2 nhóm và để trong đối tượng group > group <- cut2(bmd, g=2) > table(group) group [-3.21,-0.92) [-0.92, 2.12] 5 5 Nh thấy qua ví dụ trên, g = 2 có nghĩa là chia thành 2 nhóm (g = group). R tự động chia thành nhóm 1 gồm giá trị bmd từ -3.21 đến -0.92, và nhóm 2 từ -0.92 đến 2.12. Mỗi nhóm gồm có 5 số. Tất nhiên, chúng ta cũng có thể chia thành 3 nhóm bằng lệnh: > group <- cut2(bmd, g=3) Và với lệnh table chúng ta sẽ biết có 3 nhóm, nhóm 1 gồm 4 số, nhóm 2 và 3 mỗi nhóm có 3 số: > table(group) group [-3.21,-1.80) [-1.80, 0.21) [ 0.21, 2.12] 4 3 3 6. Sử dụng R cho tính toán đơn giản Một trong những lợi thế của R là có thể sử dụng nh một … máy tính cầm tay. Thật ra, hơn thế nữa, R có thể sử dụng cho các phép tính ma trận và lập ch ơng. Trong ch ơng này tôi chỉ trình bày một số phép tính đơn gi n mà học sinh hay sinh viên có thể sử dụng lập tức trong khi đọc những dòng chữ này. 6.1 Tính toán đơn giản 24 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Cộng hai số hay nhiều số với nhau: Cộng và trừ: > 15+2997 [1] 3012 > 15+2997-9768 [1] -6756 Nhân và chia Số lũy thừa: (25 – 5)3 > -27*12/21 [1] -15.42857 > (25 - 5)^3 [1] 8000 Căn số bậc hai: 10 Số pi (π) > sqrt(10) [1] 3.162278 > pi [1] 3.141593 > 2+3*pi [1] 11.42478 Logarit: loge Logarit: log10 > log(10) [1] 2.302585 > log10(100) [1] 2 Số mũ: e2.7689 Hàm số l ợng giác > exp(2.7689) [1] 15.94109 > cos(pi) [1] -1 > log10(2+3*pi) [1] 1.057848 Vector > x <- c(2,3,1,5,4,6,7,6,8) > x [1] 2 3 1 5 4 6 7 6 8 > exp(cos(x/10)) [1] 2.664634 2.599545 2.704736 2.405 2.511954 2.282647 2.148655 2.282647 [9] 2.007132 > sum(x) [1] 42 > x*2 [1] 4 > exp(x/10) [1] 1.221403 1.349859 1.105171 1.648 1.491825 1.822119 2.013753 1.822119 [9] 2.225541 6 2 10 8 12 14 12 16 Tính tổng bình ph ơng (sum of squares): 12 Tính tổng bình ph ơng điều chỉnh n + 22 + 32 + 42 + 52 = ? 2 (adjusted sum of squares): ∑ ( xi − x ) = ? > x <- c(1,2,3,4,5) > sum(x^2) [1] 55 i =1 > x <- c(1,2,3,4,5) > sum((x-mean(x))^2) [1] 10 Trong công thức trên mean(x) là số trung bình của vector x. Tính sai số bình ph ơng (mean square): Tính ph ơng sai (variance) và độ lệch chuẩn (standard deviation): 25 Phân tích số liệu và biểu đồ bằng R ∑( x − x ) n i =1 i 2 Nguyễn Văn Tuấn Ph ơng sai: s 2 = ∑ ( xi − x ) / ( n − 1) = ? n /n= ? 2 i =1 > x <- c(1,2,3,4,5) > sum((x-mean(x))^2)/length(x) [1] 2 Trong công thức trên, length(x) có nghĩa là tổng số phần tử (elements) trong vector x. > x <- c(1,2,3,4,5) > var(x) [1] 2.5 Độ lệch chuẩn: s2 : > sd(x) [1] 1.581139 6.2 Sử dụng R cho các phép tính ma trận Nh chúng ta biết ma trận (matrix), nói đơn gi n, gồm có dòng (row) và cột (column). Khi viết A[m, n], chúng ta hiểu rằng ma trận A có m dòng và n cột. Trong R, chúng ta cũng có thể thể hiện nh thế. Ví dụ: chúng ta muốn t o một ma trận vuông A gồm 3 dòng và 3 cột, với các phần tử (element) 1, 2, 3, 4, 5, 6, 7, 8, 9, chúng ta viết: 1 4 7   A = 2 5 8 3 6 9   Và với R: > y <- c(1,2,3,4,5,6,7,8,9) > A <- matrix(y, nrow=3) > A [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 Nh ng nếu chúng ta lệnh: > A <- matrix(y, nrow=3, byrow=TRUE) > A thì kết qu sẽ là: [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 Tức là một ma trận chuyển vị (transposed matrix). Một cách khác để t o một ma trận hoán vị là dùng t(). Ví dụ: 26 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > y <- c(1,2,3,4,5,6,7,8,9) > A <- matrix(y, nrow=3) > A [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 và B = A' có thể diễn t bằng R nh sau: > B <- t(A) > B [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 Ma trận vô hướng (scalar matrix) là một ma trận vuông (tức số dòng bằng số cột), và tất c các phần tử ngoài đ ng chéo (off-diagonal elements) là 0, và phần tử đ ng chéo là 1. Chúng ta có thể t o một ma trận nh thế bằng R nh sau: > # tạo ra mộ ma trận 3 x 3 với tất cả phần tử là 0. > A <- matrix(0, 3, 3) > # cho các phần tử đường chéo bằng 1 > diag(A) <- 1 > diag(A) [1] 1 1 1 > # bây giờ ma trận A sẽ là: > A [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1 6.2.1 Chiết phần tử từ ma trận > y <- c(1,2,3,4,5,6,7,8,9) > A <- matrix(y, nrow=3) > A [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > # cột 1 của ma trận A > A[,1] 27 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn [1] 1 4 7 > # cột 3 của ma trận A > A[3,] [1] 7 8 9 > # dòng 1 của ma trận A > A[1,] [1] 1 2 3 > # dòng 2, cột 3 của ma trận A > A[2,3] [1] 6 > # tất c các dòng của ma trận A, ngo i trừ dòng 2 > A[-2,] [,1] [,2] [,3] [1,] 1 4 7 [2,] 3 6 9 > # tất c các cột của ma trận A, ngo i trừ cột 1 > A[,-1] [,1] [,2] [1,] 4 7 [2,] 5 8 [3,] 6 9 > # xem phần tử nào cao hơn 3. > A>3 [,1] [,2] [,3] [1,] FALSE TRUE TRUE [2,] FALSE TRUE TRUE [3,] FALSE TRUE TRUE 6.2.2 Tính toán với ma trận Cộng và trừ hai ma trận. Cho hai ma trận A và B nh sau: > A <- matrix(1:12, 3, 4) > A [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12 > B <- matrix(-1:-12, 3, 4) > B [,1] [,2] [,3] [,4] [1,] -1 -4 -7 -10 28 Phân tích số liệu và biểu đồ bằng R [2,] [3,] -2 -3 -5 -6 -8 -9 Nguyễn Văn Tuấn -11 -12 Chúng ta có thể cộng A+B: > C <- A+B > C [,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 0 0 0 0 [3,] 0 0 0 0 Hay A-B: > D <- A-B > D [,1] [,2] [,3] [,4] [1,] 2 8 14 20 [2,] 4 10 16 22 [3,] 6 12 18 24 Nhân hai ma trận. Cho hai ma trận: 1 4 7   A = 2 5 8 3 6 9   và  1 2 3   B =  4 5 6 7 8 9   Chúng ta muốn tính AB, và có thể triển khai bằng R bằng cách sử dụng %*% nh sau: > > > > > y <- c(1,2,3,4,5,6,7,8,9) A <- matrix(y, nrow=3) B <- t(A) AB <- A%*%B AB [,1] [,2] [,3] [1,] 66 78 90 [2,] 78 93 108 [3,] 90 108 126 Hay tính BA, và có thể triển khai bằng R bằng cách sử dụng %*% nh sau: > BA <- B%*%A > BA [,1] [,2] [,3] [1,] 14 32 50 [2,] 32 77 122 [3,] 50 122 194 29 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Nghịch đảo ma trận và giải hệ phương trình. Ví dụ chúng ta có hệ ph ơng trình sau đây: 3x1 + 4 x2 = 4 x1 + 6 x2 = 2 Hệ ph ơng trình này có thể viết bằng kí hiệu ma trận: AX = Y, trong đó: 3 4 A= , 1 6 x  X =  1 ,  x2  và  4 Y =   2 Nghiệm của hệ ph ơng trình này là: X = A-1Y, hay trong R: > > > > A <- matrix(c(3,1,4,6), nrow=2) Y <- matrix(c(4,2), nrow=2) X <- solve(A)%*%Y X [,1] [1,] 1.1428571 [2,] 0.1428571 Chúng ta có thể kiểm tra: > 3*X[1,1]+4*X[2,1] [1] 4 Trị số eigen cũng có thể tính toán bằng function eigen nh sau: > eigen(A) $values [1] 7 2 $vectors [,1] [,2] [1,] -0.7071068 -0.9701425 [2,] -0.7071068 0.2425356 Định th c (determinant). Làm sao chúng ta xác định một ma trận có thể đ o nghịch hay không? Ma trận mà định thức bằng 0 là ma trận suy biến (singular matrix) và không thể đ o nghịch. Để kiểm tra định thức, R dùng lệnh det(): > E <- matrix((1:9), 3, 3) > E [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 30 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > det(E) [1] 0 Nh ng ma trận F sau đây thì có thể đ o nghịch: > F <- matrix((1:9)^2, 3, 3) > F [,1] [,2] [,3] [1,] 1 16 49 [2,] 4 25 64 [3,] 9 36 81 > det(F) [1] -216 Và nghịch đ o của ma trận F (F-1) có thể tính bằng function solve() nh sau: > solve(F) [,1] [,2] [,3] [1,] 1.291667 -2.166667 0.9305556 [2,] -1.166667 1.666667 -0.6111111 [3,] 0.375000 -0.500000 0.1805556 Ngoài những phép tính đơn gi n này, R còn có thể sử dụng cho các phép tính phức t p khác. Một lợi thế đáng kể của R là phần mềm cung cấp cho ng i sử dụng tự do t o ra những phép tính phù hợp cho từng vấn đề cụ thể. R có một package Matrix chuyên thiết kế cho tính toán ma trận. B n đọc có thể t i package xuống, cài vào máy, và sử dụng, nếu cần. Địa chỉ để t i là: http://cran.au.r-project.org/bin/windows/contrib/rrelease/Matrix_0.995-8.zip cùng với tài liệu chỉ dẫn cách sử dụng (dài kho ng 80 trang): http://cran.au.r-project.org/doc/packages/Matrix.pdf. 7. Sử dụng R cho tính toán xác suất 7.1 Phép hoán vị (permutation) Chúng ta biết rằng 3! = 3.2.1 = 6, và 0!=1. Nói chung, công thức tính hoán vị cho một số n là: n ! = n ( n − 1)( n − 2 )( n − 3) × ... × 1 . Trong R cách tính này rất đơn gi n với lệnh prod() nh sau: • Tìm 3! > prod(3:1) [1] 6 • Tìm 10! > prod(10:1) [1] 3628800 31 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn • Tìm 10.9.8.7.6.5.4 > prod(10:4) [1] 604800 • Tìm (10.9.8.7.6.5.4) / (40.39.38.37.36) > prod(10:4) / prod(40:36) [1] 0.007659481 7.2 Tổ hợp (combination) Số lần chọn k ng n n! i từ n phần tử là:   = . Công thức này cũng có khi viết là  k  k !( n − k ) ! n Ckn thay vì   . Với R, phép tính này rất đơn gi n bằng hàm choose(n, k). Sau k  đây là vài ví dụ minh họa: 5 Tìm    2 > choose(5, 2) [1] 10 • • Tìm xác suất cặp A và B trong số 5 ng > 1/choose(5, 2) [1] 0.1 i đ ợc đắc cử vào hai chức vụ: 7.3 Biến số ngẫu nhiên và hàm phân phối Khi nói đến “phân phối” (hay distribution) là đề cập đến các giá trị mà biến số có thể có. Các hàm phân phối (distribution function) là hàm nhằm mô t các biến số đó một cách có hệ thống. “Có hệ thống” đây có nghĩa là theo mộ mô hình toán học cụ thể với những thông số cho tr ớc. Trong xác suất thống kê có khá nhiều hàm phân phối, và đây chúng ta sẽ xem xét qua một số hàm quan trọng nhất và thông dụng nhất: đó là phân phối nhị phân, phân phối Poisson, và phân phối chuẩn. Trong mỗi luật phân phối, có 4 lo i hàm quan trọng mà chúng ta cần biết: • • • • hàm mật độ xác suất (probability density distribution); hàm phân phối tích lũy (cumulative probability distribution); hàm định bậc (quantile); và hàm mô phỏng (simulation). R có những hàm sẵn trên có thể ứng dụng cho tính toán xác suất. Tên mỗi hàm đ ợc gọi bằng một tiếp đầu ngữ để chỉ lo i hàm phân phối, và viết tắt tên của hàm đó. Các tiếp đầu ngữ là d (chỉ distribution hay xác suất), p (chỉ cumulative probability, xác suất tích lũy), q (chỉ định bậc hay quantile), và r (chỉ random hay số ngẫu nhiên). Các 32 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn tên viết tắt là norm (normal, phân phối chuẩn), binom (binomial , phân phối nhị phân), pois (Poisson, phân phối Poisson), v.v… B ng sau đây tóm tắt các hàm và thông số cho từng hàm: Hàm phân phối Mật độ Tích lũy Định bậc Mô phỏng Chuẩn dnorm(x, mean, sd) dbinom(k, n, p) pnorm(q, mean, sd) qnorm(p, mean, sd) rnorm(n, mean, sd) pbinom(q, n, p) qbinom (p, n, p) rbinom(k, n, prob) dpois(k, lambda) ppois(q, lambda) qpois(p, lambda) rpois(n, lambda) dunif(x, min, max) dnbinom(x, k, p) punif(q, min, max) qunif(p, min, max) runif(n, min, max) pnbinom(q, k, p) qnbinom (p,k,prob) rbinom(n, n, prob) dbeta(x, shape1, shape2) dgamma(x, shape, rate, scale) dgeom(x, p) pbeta(q, shape1, shape2) gamma(q, shape, rate, scale) pgeom(q, p) qbeta(p, shape1, shape2) qgamma(p, shape, rate, scale) qgeom(p, prob) rbeta(n, shape1, shape2) rgamma(n, shape, rate, scale) rgeom(n, prob) dexp(x, rate) pexp(q, rate) qexp(p, rate) rexp(n, rate) dnorm(x, mean, sd) dcauchy(x, location, scale) df(x, df1, df2) pnorm(q, mean, sd) qnorm(p, mean, sd) rnorm(n, mean, sd) pcauchy(q, location, scale) pf(q, df1, df2) qcauchy(p, location, scale) qf(p, df1, df2) rcauchy(n, location, scale) rf(n, df1, df2) Nhị phân Poisson Uniform Negative binomial Beta Gamma Geometric Exponential Weibull Cauchy F dt(x, df) pt(q, df) qt(p, df) rt(n, df) T dchisq(x, df) pchi(q, df) qchisq(p, df) rchisq(n, df) Chi-squared Chú thích: Trong b ng trên, df = degrees of freedome (bậc tự do); prob = probability (xác suất); n = sample size (số l ợng mẫu). Các thông số khác có thể tham kh o thêm cho từng luật phân phối. Riêng các luật phân phối F, t, Chi-squared còn có một thông số khác nữa là non-centrality parameter (ncp) đ ợc cho số 0. Tuy nhiên ng i sử dụng có thể cho một thông số khác thích hợp, nếu cần. 7.3.1 Hàm phân phối nhị phân (Binomial distribution) Nh tên gọi, hàm phân phối nhị phân chỉ có hai giá trị: nam / nữ, sống / chết, có / không, v.v… Hàm nhị phân đ ợc phát biểu bằng định lí nh sau: Nếu một thử nghiệm đ ợc tiến hành n lần, mỗi lần cho ra kết qu hoặc là thành công hoặc là thất b i, và gồm xác suất thành công đ ợc biết tr ớc là p, thì xác suất có k lần thử nghiệm thành công là: n−k P ( k | n, p ) = Ckn p k (1 − p ) , trong đó k = 0, 1, 2, . . . , n. Trong R, có hàm dbinom(k, n, p) có thể giúp chúng ta tính công thức P ( k | n, p ) = Ckn p k (1 − p ) chóng. Trong tr n−k một cách nhanh ng hợp trên, chúng ta chỉ cần đơn gi n lệnh: > dbinom(2, 3, 0.60) [1] 0.432 Ví dụ 2: Hàm nhị phân tích lũy (Cumulative Binomial probability distribution). Xác suất thuốc chống loãng x ơng có hiệu nghiệm là kho ng 70% (tức là p = 0.70). Nếu chúng ta điều trị 10 bệnh nhân, xác suất có tối thiểu 8 bệnh nhân với kết qu tích cực là bao nhiêu? Nói cách khác, nếu gọi X là số bệnh nhân đ ợc điều trị thành công, chúng ta cần tìm P(X ≥ 8) = ? Để tr l i câu hỏi này, chúng ta sử dụng hàm 33 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn pbinom(k, n, p). Xin nhắc l i rằng hàm pbinom(k, n, p)cho chúng ta P(X ≤ k). Do đó, P(X ≥ 8) = 1 – P(X ≤ 7). Thành ra, đáp số bằng R cho câu hỏi là: > 1-pbinom(7, 10, 0.70) [1] 0.3827828 Ví dụ 3: Mô phỏng hàm nhị phân: Biết rằng trong một quần thể dân số có kho ng 20% ng i mắc bệnh cao huyết áp; nếu chúng ta tiến hành chọn mẫu 1000 lần, mỗi lần chọn 20 ng i trong quần thể đó một cách ngẫu nhiên, sự phân phối số bệnh nhân cao huyết áp sẽ nh thế nào? Để tr l i câu hỏi này, chúng ta có thể ứng dụng hàm rbinom (n, k, p) trong R với những thông số nh sau: > b <- rbinom(1000, 20, 0.20) Trong lệnh trên, kết qu mô phỏng đ ợc t m th i chứa trong đối t ợng tên là b. Để biết b có gì, chúng ta đếm bằng lệnh table: > table(b) b 0 1 2 3 4 5 6 6 45 147 192 229 169 105 7 68 8 23 9 13 10 3 Dòng số liệu thứ nhất (0, 5, 6, …, 10) là số bệnh nhân mắc bệnh cao huyết áp trong số 20 ng i mà chúng ta chọn. Dòng số liệu thứ hai cho chúng ta biết số lần chọn mẫu trong 1000 lần x y ra. Do đó, có 6 mẫu không có bệnh nhân cao huyết áp nào, 45 mẫu với chỉ 1 bệnh nhân cao huyết áp, v.v… Có lẽ cách để hiểu là vẽ đồ thị các tần số trên bằng lệnh hist nh sau: > hist(b, main="Number of hypertensive patients") 0 50 Frequency 100 150 200 Number of hypertensive patients 0 2 4 6 8 10 b 34 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Biểu đồ 1. Phân phối số bệnh nhân cao huyết áp trong số 20 ng i đ ợc chọn ngẫu nhiên trong một quần thề gồm 20% bệnh nhân cao huyết áp, và chọn mẫu đ ợc lặp l i 1000 lần. Qua biểu đồ trên, chúng ta thấy xác suất có 4 bệnh nhân cao huyết áp (trong mỗi lần chọn mẫu 20 ng i) là cao nhất (22.9%). Điều này cũng có thể hiểu đ ợc, b i vì tỉ lệ cao huyết áp là 20%, cho nên chúng ta kì vọng rằng trung bình 4 ng i trong số 20 ng i đ ợc chọn ph i là cao huyết áp. Tuy nhiên, điều quan trọng mà biểu đồ trên thể hiện là có khi chúng ta quan sát đến 10 bệnh nhân cao huyết áp dù xác suất cho mẫu này rất thấp (chỉ 3/1000). 7.3.2 Hàm phân phối Poisson (Poisson distribution) Hàm phân phối Poisson, nói chung, rất giống với hàm nhị phân, ngo i trừ thông số p th ng rất nhỏ và n th ng rất lớn. Vì thế, hàm Poisson th ng đ ợc sử dụng để mô t các biến số rất hiếm x y ra (nh số ng i mắc ung th trong một dân số chẳng h n). Hàm Poisson còn đ ợc ứng dụng khá nhiều và thành công trong các nghiên cứu kĩ thuật và thị tr ng nh số l ợng khách hàng đến một nhà hàng mỗi gi . Ví dụ 4: Hàm mật độ Poisson (Poisson density probability function). Qua theo dõi nhiều tháng, ng i ta biết đ ợc tỉ lệ đánh sai chính t của một th kí đánh máy. Tính trung bình cứ kho ng 2.000 chữ thì th kí đánh sai 1 chữ. Hỏi xác suất mà th kí đánh sai chính t 2 chữ, hơn 2 chữ là bao nhiêu? Vì tần số khá thấp, chúng ta có thể gi định rằng biến số “sai chính t ” (t m đặt tên là biến số X) là một hàm ngẫu nhiên theo luật phân phối Poisson. đây, chúng ta có tỉ lệ sai chính t trung bình là 1(λ = 1). Luật phân phối Poisson phát biểu rằng xác suất mà X = k, với điều kiện tỉ lệ trung bình λ, : P( X = k | λ) = e− λ λ k k! e−212 Do đó, đáp số cho câu hỏi trên là: P ( X = 2 | λ = 1) = = 0.1839 . Đáp số này có thể 2! tính bằng R một cách nhanh chóng hơn bằng hàm dpois nh sau: > dpois(2, 1) [1] 0.1839397 Chúng ta cũng có thể tính xác suất sai 1 chữ, và xác suất không sai chữ nào: > dpois(1, 1) [1] 0.3678794 > dpois(0, 1) 35 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn [1] 0.3678794 Chú ý trong hàm trên, chúng ta chỉ đơn gi n cung cấp thông số k = 2 và (λ = 1. Trên đây là xác suất mà th kí đánh sai chính t đúng 2 chữ. Nh ng xác suất mà th kí đánh sai chính t hơn 2 chữ (tức 3, 4, 5, … chữ) có thể ớc tính bằng: P ( X > 2 ) = P ( X = 3) + P ( X = 4 ) + P( X = 5) + ... = 1 − P ( X ≤ 2) = 1 – 0.3678 – 0.3678 – 0.1839 = 0.08 Bằng R, chúng ta có thể tính nh sau: # P(X ≤ 2) > ppois(2, 1) [1] 0.9196986 # 1-P(X ≤ 2) > 1-ppois(2, 1) [1] 0.0803014 7.3.3 Hàm phân phối chuẩn (Normal distribution) Hai luật phân phối mà chúng ta vừa xem xét trên đây thuộc vào nhóm phân phối áp dụng cho các biến số phi liên tục (discrete distributions), mà trong đó biến số có những giá trị theo bậc thứ hay thể lo i. Đối với các biến số liên tục, có vài luật phân phối thích hợp khác, mà quan trọng nhất là phân phối chuẩn. Phân phối chuẩn là nền t ng quan trọng nhất của phân tích thống kê. Có thể nói không ngoa rằng hầu hết lí thuyết thống kê đ ợc xây dựng trên nền t ng của phân phối chuẩn. Hàm mật độ phân phối chuẩn có hai thông số: trung bình µ và ph ơng sai σ2 (hay độ lệch chuẩn σ). Gọi X là một biến số (nh chiều cao chẳng h n), hàm mật độ phân phối chuẩn phát biểu rằng xác suất mà X = x là:  ( x − µ )2  1 2 P ( X = x | µ ,σ ) = f ( x ) = exp  −  2σ 2  2πσ   Ví dụ 5: Hàm mật độ phân phối chuẩn (Normal density probability function). Chiều cao trung bình hiện nay phụ nữ Việt Nam là 156 cm, với độ lệch chuẩn là 4.6 cm. Cũng biết rằng chiều cao này tuân theo luật phân phối chuẩn. Với hai thông số µ=156, σ=4.6, chúng ta có thể xây dựng một hàm phân phối chiều cao cho toàn bộ quần thể phụ nữ Việt Nam, và hàm này có hình d ng nh sau: 36 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn f(height) 0.00 0.02 0.04 0.06 0.08 Probability distribution of height in Vietnamese women 130 140 150 160 170 180 190 200 Height Biểu đồ 2. Phân phối chiều cao phụ nữ Việt Nam với trung bình 156 cm và độ lệch chuẩn 4.6 cm. Trụng hoành là chiều cao và trục tung là xác suất cho mỗi chiều cao. Biểu đồ trên đ ợc vẽ bằng hai lệnh sau đây. Lệnh đầu tiên nhằm t o ra một biến số height có giá trị 130, 131, 132, …, 200 cm. Lệnh thứ hai là vẽ biểu đồ với điều kiện trung bình là 156 cm và độ lệch chuẩn là 4.6 cm. > height <- seq(130, 200, 1) > plot(height, dnorm(height, 156, 4.6), type="l", ylab=”f(height)”, xlab=”Height”, main="Probability distribution of height in Vietnamese women") Với hai thông số trên (và biểu đồ), chúng ta có thể ớc tính xác suất cho bất cứ chiều cao nào. Chẳng h n nh xác suất một phụ nữ Việt Nam có chiều cao 160 cm là:  (160 − 156 )2  1 P(X = 160 | µ=156, σ=4.6) = exp  −  2 4.6 2 × 3.1416 2 ( 4.6 )   = 0.0594 Hàm dnorm(x, mean, sd)trong R có thể tính toán xác suất này cho chúng ta một cách gọn nhẹ: > dnorm(160, mean=156, sd=4.6) [1] 0.05942343 37 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Hàm xác suất chuẩn tích lũy (cumulative normal probability function). Vì chiều cao là một biến số liên tục, trong thực tế chúng ta ít khi nào muốn tìm xác suất cho một giá trị cụ thể x, mà th ng tìm xác suất cho một kho ng giá trị a đến b. Chẳng h n nh chúng ta muốn biết xác suất chiều cao từ 150 đến 160 cm (tức là P(160 ≤ X ≤ 150), hay xác suất chiều cao thấp hơn 145 cm, tức P(X < 145). Để tìm đáp số các câu hỏi nh thế, chúng ta cần đến hàm xác suất chuẩn tích lũy, đ ợc định nghĩa nh sau: P(a ≤ X ≤ b) = ∫ f ( x ) dx b a Thành ra, P(160 ≤ X ≤ 150) chính là diện tích tính từ trục hoành = 150 đến 160 của biểu đồ 2. Trong R có hàm pnorm(x, mean, sd) dùng để tính xác suất tích lũy cho một phân phối chuẩn rất có ích. pnorm (a, mean, sd) = ∫ a −∞ f ( x ) dx = P(X ≤ a | mean, sd) Chẳng h n nh xác suất chiều cao phụ nữ Việt Nam bằng hoặc thấp hơn 150 cm là 9.6%: > pnorm(150, 156, 4.6) [1] 0.0960575 Hay xác suất chiều cao phụ nữ Việt Nam bằng hoặc cao hơn 165 cm là: > 1-pnorm(164, 156, 4.6) [1] 0.04100591 Nói cách khác, chỉ có kho ng 4.1% phụ nữ Việt Nam có chiều cao bằng hay cao hơn 165 cm. Ví dụ 6: ng dụng luật phân phối chuẩn: Trong một quần thể, chúng ta biết rằng áp suất máu trung bình là 100 mmHg và độ lệch chuẩn là 13 mmHg, hỏi: có bao nhiêu ngừơi trong quần thể này có áp suất máu bằng hoặc cao hơn 120 mmHg? Câu tr l i bằng R là: > 1-pnorm(120, mean=100, sd=13) [1] 0.0619679 Tức kho ng 6.2% ng mmHg. i trong quần thể này có áp suất máu bằng hoặc cao hơn 120 7.3.4 Hàm phân phối chuẩn chuẩn hóa (Standardized Normal distribution) 38 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Một biến X tuân theo luật phân phối chuẩn với trung bình bình µ và ph ơng sai σ2 ng đ ợc viết tắt là: X ~ N(µ , σ2) th đây µ và σ2 tùy thuộc vào đơn vị đo l ng của biến số. Chẳng h n nh chiều cao đ ợc tính bằng cm (hay m), huyết áp đ ợc đo bằng mmHg, tuổi đ ợc đo bằng năm, v.v… cho nên đôi khi mô t một biến số bằng đơn vị gốc rất khó so sánh. Một cách đơn gi n hơn là chuẩn hóa (standardized) X sao cho số trung bình là 0 và ph ơng sai là 1. Sau vài thao tác số học, có thể chứng minh dễ dàng rằng, cách biến đổi X để đáp ứng điều kiện trên là: X −µ Z= σ Nói theo ngôn ngữ toán: nếu X ~ N(µ , σ2), thì (X – µ)/σ2 ~ N(0, 1). Nh vậy qua công thức trên, Z thực chất là độ khác biệt giữa một số và trung bình tính bằng số độ lệch chuẩn. Nếu Z = 0, chúng ta biết rằng X bằng số trung bình µ. Nếu Z = -1, chúng ta biết rằng X thấp hơn µ đúng 1 độ lệch chuẩn. T ơng tự, Z = 2.5, chúng ta biết rằng X cao hơn µ đúng 2.5 độ lệch chuẩn. v.v… Biểu đồ phân phối chiều cao của phụ nữ Việt Nam có thể mô t bằng một đơn vị mới, đó là chỉ số z nh sau: 0.2 0.0 0.1 f(z) 0.3 0.4 Probability distribution of height in Vietnamese women -4 -2 0 2 4 z Biểu đồ 3. Phân phối chuẩn hóa chiều cao nữ Việt Nam. phụ Biểu đồ trên đ ợc vẽ bằng hai lệnh sau đây: 39 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn > height <- seq(-4, 4, 0.1) > plot(height, dnorm(height, 0, 1), type="l", ylab=”f(z)”, xlab=”z”, main="Probability distribution of height in Vietnamese women") Với phân phối chuẩn chuẩn hoá, chúng ta có một tiện lợi là có thể dùng nó để mô t và so sánh mật độ phân phối của bất cứ biến nào, vì tất c đều đ ợc chuyển sang chỉ số z. Trong biểu đồ trên, trục tung là xác suất z và trục hoành là biến số z. Chúng ta có thể tính toán xác suất z nhỏ hơn một hằng số (constant) nào đó dê dàng bằng R. Ví dụ, chúng ta muốn tìm P(z ≤ -1.96) = ? cho một phân phối mà trung bình là 0 và độ lệch chuẩn là 1. > pnorm(-1.96, mean=0, sd=1) [1] 0.02499790 Hay P(z ≤ 1.96) = ? > pnorm(1.96, mean=0, sd=1) [1] 0.9750021 Do đó, P(-1.96 < z < 1.96) chính là: > pnorm(1.96) - pnorm(-1.96) [1] 0.9500042 Nói cách khác, xác suất 95% là z nằm giữa -1.96 và 1.96. (Chú ý trong lệnh trên tôi không cung cấp mean=0, sd=1, b i vì trong thực tế, pnorm giá trị mặc định (default value) của thông số mean là 0 và sd là 1). Ví dụ 5 (tiếp tục). Xin nhắc l i để tiện việc theo dõi, chiều cao trung bình phụ nữ Việt Nam là 156 cm và độ lệch chuẩn là 4.6 cm. Do đó, một phụ nữ có chiều cao 170 cm cũng có nghĩa là z = (170 – 156) / 4.6 = 3.04 độ lệch chuẩn, và ti lệ các phụ nữ Việt Nam có chiều cao cao hơn 170 cm là rất thấp, chỉ kho ng 0.1%. > 1-pnorm(3.04) [1] 0.001182891 Tìm định lượng (quantile) c a một phân phối chuẩn. Đôi khi chúng ta cần làm một tính toán đ o ng ợc. Chẳng h n nh chúng ta muốn biết: nếu xác suất Z nhỏ hơn một hằng số z nào đó cho tr ớc bằng p, thì z là bao nhiêu? Diễn t theo kí hiệu xác suất, chúng ta muốn tìm z trong nếu: P(Z < z) = p Để tr l i câu hỏi này, chúng ta sử dụng hàm qnorm(p, mean=, sd=). 40 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ví dụ 7: Biết rằng Z ~ N(0, 1) và nếu P(Z < z) = 0.95, chúng ta muốn tìm z. > qnorm(0.95, mean=0, sd=1) [1] 1.644854 Hay P(Z < z) = 0.975 cho phân phối chuẩn với trung bình 0 và độ lệch chuẩn 1: > qnorm(0.975, mean=0, sd=1) [1] 1.959964 7.4 Chọn mẫu ngẫu nhiên (random sampling) Trong xác suất và thống kê, lấy mẫu ngẫu nhiên rất quan trọng, vì nó đ m b o tính hợp lí của các ph ơng pháp phân tích và suy luận thống kê. Với R, chúng ta có thể lấy mẫu một mẫu ngẫu nhiên bằng cách sử dụng hàm sample. Ví dụ 8: Chúng ta có một quần thể gồm 40 ng i (mã số 1, 2, 3, …, 40). Nếu chúng ta muốn chọn 5 đối t ợng quần thể đó, ai sẽ là ng i đ ợc chọn? Chúng ta có thể dùng lệnh sample() để tr l i câu hỏi đó nh sau: > sample(1:40, 5) [1] 32 26 6 18 9 Kết qu trên cho biết đối t ợng 32, 26, 8, 18 và 9 đ ợc chọn. Mỗi lần ra lệnh này, R sẽ chọn một mẫu khác, chứ không hoàn toàn giống nh mẫu trên. Ví dụ: > sample(1:40, 5) [1] 5 22 35 19 4 > sample(1:40, 5) [1] 24 26 12 6 22 > sample(1:40, 5) [1] 22 38 11 6 18 v.v… Trên đây là lệnh để chúng ta chọn mẫu ngẫu nhiên mà không thay thế (random sampling without replacement), tức là mỗi lần chọn mẫu, chúng ta không bỏ l i các mẫu đã chọn vào quần thể. Nh ng nếu chúng ta muốn chọn mẫu thay thế (tức mỗi lần chọn ra một số đối t ợng, chúng ta bỏ vào l i trong quần thể để chọn tiếp lần sau). Ví dụ, chúng ta muốn chọn 10 ng i từ một quần thể 50 ng i, bằng cách lấy mẫu với thay thế (random sampling with replacement), chúng ta chỉ cần thêm tham số replace = TRUE: > sample(1:50, 10, replace=T) 41 Phân tích số liệu và biểu đồ bằng R [1] 31 44 6 Nguyễn Văn Tuấn 8 47 50 10 16 29 23 Hay ném một đồng xu 10 lần; mỗi lần, dĩ nhiên đồng xu có 2 kết qu H và T; và kết qu 10 lần có thể là: > sample(c("H", "T"), 10, replace=T) [1] "H" "T" "H" "H" "H" "T" "H" "H" "T" "T" Cũng có thể t ng t ợng chúng ta có 5 qu banh màu xanh (X) và 5 qu banh màu đỏ (D) trong một bao. Nếu chúng ta chọn 1 qu banh, ghi nhận màu, rồi để l i vào bao; rồi l i chọn 1 qu banh khác, ghi nhận màu, và bỏ vào bao l i. Cứ nh thế, chúng ta chọn 20 lần, kết qu có thể là: > sample(c("X", "D"), 20, replace=T) [1] "X" "D" "D" "D" "D" "D" "X" "X" "X" "X" "X" "D" "X" "X" "D" "X" "X" "X" "X" [20] "D" Ngoài ra, chúng ta còn có thể lấy mẫu với một xác suất cho tr ớc. Trong hàm sau đây, chúng ta chọn 10 đối t ợng từ dãy số 1 đến 5, nh ng xác suất không bằng nhau: > sample(5, 10, prob=c(0.3, 0.4, 0.1, 0.1, 0.1), replace=T) [1] 3 1 3 2 2 2 2 2 5 1 Đối t ợng 1 đ ợc chọn 2 lần, đối t ợng 2 đ ợc chọn 5 lần, đối t ợng 3 đ ợc chọn 2 lần, v.v… Tuy không hoàn toàn phù hợp với xác suất 0.3, 0.4, 0.1 nh cung cấp vì số mẫu còn nhỏ, nh ng cũng không quá xa với kì vọng. 8. Biểu đồ Trong ngôn ngữ R có rất nhiều cách để thiết kế một biểu đồ gọn và đẹp. Phần lớn những hàm để thiết kế biểu đồ có sẵn trong R, nh ng một số lo i biểu đồ tinh vi và phức t p khác có thể thiết kế bằng các package chuyên dụng nh lattice hay trellis có thể t i từ website của R. Trong ch ơng này tôi sẽ chỉ cách vẽ các biểu đồ thông dụng bằng cách sử dụng các hàm phổ biến trong R. 8.1 Số liệu cho phân tích biểu đồ Sau khi đã biết qua môi tr ng và những lựa chọn để thiết kế một biểu đồ, bây gi chúng ta có thể sử dụng một số hàm thông dụng để vẽ các biểu đồ cho số liệu. Theo tôi, biểu đồ có thể chia thành 2 lo i chính: biểu đồ dùng để mô t một biến số và biểu đồ về mối liên hệ giữa hai hay nhiều biến số. Tất nhiên, biến số có thể là liên tục hay không liên tục, cho nên, trong thực tế, chúng ta có 4 lo i biểu đồ. Trong phần sau đây, tôi sẽ điểm qua các lo i biểu đồ, từ đơn gi n đến phức t p. Có lẽ cách tốt nhất để tìm hiểu cách vẽ đồ thị bằng R là bằng một dữ liệu thực tế. Tôi sẽ quay l i ví dụ 2 (phần 4.2). Trong ví dụ đó, chúng ta có dữ liệu gồm 8 cột (hay 42 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn biến số): id, sex, age, bmi, hdl, ldl, tc, và tg. (Chú ý, id là mã số của 50 đối t ợng nghiên cứu; sex là giới tính (nam hay nữ); age là độ tuổi; bmi là tỉ số trọng l ơng; hdl là high density cholesterol; ldl là low density cholesterol; tc là tổng số - total – cholesterol; và tg triglycerides). Dữ liệu đ ợc chứa trong directory directory c:\works\insulin d ới tên chol.txt. Tr ớc khi vẽ đồ thị, chúng ta bắt đầu bằng cách nhập dữ liệu này vào R. > setwd(“c:/works/stats”) > cong <- read.table(“chol.txt”, header=TRUE, na.strings=”.”) > attach(cong) Hay để tiện việc theo dõi tôi sẽ nhập các dữ liệu đó bằng các lệnh sau đây: sex <- c(“Nam”, “Nu”, “Nu”,“Nam”,“Nam”, “Nu”,“Nam”,“Nam”,“Nam”, “Nu”, “Nu”,“Nam”, “Nu”,“Nam”,“Nam”, “Nu”, “Nu”, “Nu”, “Nu”, “Nu”, “Nu”, “Nu”, “Nu”, “Nu”,“Nam”,“Nam”, “Nu”,“Nam”, “Nu”, “Nu”, “Nu”,“Nam”,“Nam”, “Nu”, “Nu”,“Nam”, “Nu”,“Nam”, “Nu”, “Nu”, “Nam”, “Nu”,“Nam”,“Nam”,“Nam”, “Nu”,“Nam”,“Nam”, “Nu”, “Nu”) age <- c(57, 63, 61, 60, 51, 64, 51, 45, 50, 58, bmi <- c( 17, 20, 22, 24, 60, 60, 70, 60, 60, 18, 21, 22, 24, 65, 42, 51, 55, 45, 18, 21, 22, 24, 47, 64, 63, 74, 63, 18, 21, 22, 25, 65, 49, 54, 48, 52, 76, 44, 57, 46, 64, 61, 45, 70, 49, 45, 59, 80, 47, 69, 64, 57, 48, 60, 72, 62) 18, 18, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 25) hdl <- c(5.000,4.380,3.360,5.920,6.250,4.150,0.737,7.170,6.942,5.000, 4.217,4.823,3.750,1.904,6.900,0.633,5.530,6.625,5.960,3.800, 5.375,3.360,5.000,2.608,4.130,5.000,6.235,3.600,5.625,5.360, 6.580,7.545,6.440,6.170,5.270,3.220,5.400,6.300,9.110,7.750, 6.200,7.050,6.300,5.450,5.000,3.360,7.170,7.880,7.360,7.750) ldl <- c(2.0, 5.0, 3.1, 4.4, 3.0, 3.0, 1.3, 3.0, 4.3, 4.1, 3.0, 1.2, 1.7, 2.3, 4.4, 4.0, 0.7, 2.0, 6.0, 2.8, 2.1, 4.0, 2.1, 3.0, 3.0, 3.0, 4.1, 4.0, 3.0, 2.0, 3.0, 4.3, 4.1, 2.6, 1.0, 3.0, 4.0, 4.0, 4.4, 4.0, 3.0, 4.3, 4.2, 4.3, 4.6, tc <-c (4.0, 6.2, 4.3, 5.6, 6.2, 3.5, 4.1, 4.8, 8.3, 6.7, 4.7, 3.0, 4.0, 5.8, 6.3, 7.7, 4.0, 3.0, 7.6, 6.0, 5.0, 6.9, 3.1, 5.8, 4.0, 4.2, 5.7, 5.3, 3.1, 3.7, 5.9, 5.7, 5.3, 5.4, 6.1, 6.1, 5.3, 5.4, 6.3, 6.7, 5.9, 7.1, 4.5, 8.2, 8.1, 4.0, 3.8, 5.9, 6.2, 6.2) tg <- c(1.1, 1.7, 2.2, 3.3, 2.4, 2.1, 1.0, 2.7, 3.0, 3.3, 0.8, 1.6, 1.1, 1.0, 2.0, 1.1, 1.1, 0.7, 1.4, 2.6, 2.1, 1.5, 1.0, 2.5, 1.8, 1.5, 1.0, 1.7, 0.7, 1.2, 2.6, 2.7, 2.9, 2.4, 1.9, 1.5, 3.9, 2.5, 2.4, 3.3, 5.4, 3.0, 6.2, 1.4, 4.0, 1.9, 3.1, 1.3, 2.7, 2.5) cong <- data.frame(sex, age, bmi, hdl, ldl, tc, tg) 43 2.0, 4.0, 4.2, 4.0, 4.0) Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 8.2 Biểu đồ cho một biến số rời rạc (discrete variable): barplot Biến sex trong dữ liệu trên có hai giá trị (nam và nu), tức là một biến không liên tục. Chúng ta muốn biết tần số của giới tính (bao nhiêu nam và bao nhiêu nữ) và vẽ một biểu đồ đơn gi n. Để thực hiện ý định này, tr ớc hết, chúng ta cần dùng hàm table để biết tần số: > sex.freq <- table(sex) > sex.freq sex Nam Nu 22 28 Có 22 nam và 28 nữa trong nghiên cứu. Sau đó dùng hàm barplot để thể hiện tần số này nh sau: > barplot(sex.freq, main=”Frequency of males and females”) Biểu trên cũng có thể có đ ợc bằng một lệnh đơn gi n hơn (Biểu đồ 8a): > barplot(table(sex), main=”Frequency of males and females”) Frequency of males and females 0 5 Nam 10 15 20 Nu 25 Frequency of males and females Nam Nu Biểu đồ 8a. Tần số giới tính thể hiện bằng cột số. 0 5 10 15 20 25 Biểu đồ 8b. Tần số giới tính thể hiện bằng dòng số. Thay vì thể hiện tần số nam và nữ bằng 2 cột, chúng ta có thể thể hiện bằng hai dòng bằng thông số horiz = TRUE, nh sau (xem kết qu trong Biểu đồ 6b): > barplot(sex.freq, horiz = TRUE, col = rainbow(length(sex.freq)), main=”Frequency of males and females”) 44 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 8.3 Biểu đồ cho hai biến số rời rạc (discrete variable): barplot Age là một biến số liên tục. Chúng ta có thể chia bệnh nhân thành nhiều nhóm dựa vào độ tuổi. Hàm cut có chức năng “cắt” một biến liên tục thành nhiều nhóm r i r c. Chẳng h n nh : > ageg <- cut(age, 3) > table(ageg) ageg (42,54.7] (54.7,67.3] 19 24 (67.3,80] 7 Có hiệu qu chia biến age thành 3 nhóm. Tần số của ba nhóm này là: 42 tuổi đến 54.7 tuổi thành nhóm 1, 54.7 đến 67.3 thành nhóm 2, và 67.3 đến 80 tuổi thành nhóm 3. Nhóm 1 có 19 bệnh nhân, nhóm 2 và 3 có 24 và 7 bệnh nhân. Bây gi chúng ta muốn biết có bao nhiêu bệnh nhân trong từng độ tuổi và từng giới tính bằng lệnh table: > age.sex <- table(sex, ageg) > age.sex ageg sex (42,54.7] (54.7,67.3] (67.3,80] Nam 10 10 2 Nu 9 14 5 Kết qu trên cho thấy chúng ta có 10 bệnh nhân nam và 9 nữ trong nhóm tuổi thứ nhất, 10 nam và 14 nữa trong nhóm tuổi thứ hai, v.v… Để thể hiện tần số của hai biến này, chúng ta vẫn dùng barplot: > barplot(age.sex, main=”Number of males and females in each age group”) 45 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 0 0 2 5 4 10 6 8 15 10 20 12 14 Number of males and females in each age group (42,54.7] (54.7,67.3] (42,54.7] (67.3,80] (54.7,67.3] (67.3,80] Age group Biểu đồ 9b. Tần số giới tính và nhóm tuổi thể hiện bằng hai dòng số. Biểu đồ 9a. Tần số giới tính và nhóm tuổi thể hiện bằng cột số. Trong Biểu đồ 9a, mỗi cột là cho một độ tuổi, và phần đậm của cột là nữ, và phần màu nh t là tần số của nam giới. Thay vì thể hiện tần số nam nữ trong một cột, chúng ta cũng có thể thể hiện bằng 2 cột với beside=T nh sau (Biểu đồ 9b): barplot(age.sex, beside=TRUE, xlab="Age group") 8.4 Biểu đồ hình tròn Tần số một biến r i r c cũng có thể thể hiện bằng biểu đồ hình tròn. Ví dụ sau đây vẽ biểu đồ tần số của độ tuổi. Biểu đồ 10a là 3 nhóm độ tuổi, và Biểu đồ 10b là biểu đồ tần số cho 5 nhóm tuổi: > pie(table(ageg)) pie(table(cut(age,5))) 46 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn (42,54.7] (49.6,57.2] (42,49.6] (72.4,80] (67.3,80] (64.8,72.4] (54.7,67.3] (57.2,64.8] Biểu đồ 10a. Tần số cho 3 nhóm tuổi Biểu đồ 10b. Tần số cho 5 nhóm tuổi 8.5 Biểu đồ cho một biến số liên tục: stripchart và hist 8.5.1 Stripchart Biểu đồ strip cho chúng ta thấy tính liên tục của một biến số. Chẳng h n nh chúng ta muốn tìm hiểu tính liên tục của triglyceride (tg), hàm stripchart() sẽ giúp trong mục tiêu này: > stripchart(tg, main=”Strip chart for triglycerides”, xlab=”mg/L”) Strip chart for triglycerides 1 2 3 4 5 6 mg/L 47 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Chúng ta thấy biến số tg có sự bất liên tục, nhất là các đối t ợng có tg cao. Trong khi phần lớn đối t ợng có độ tg thấp hơn 5, thì có 2 đối t ợng với tg rất cao (>5). 8.5.2 Histogram Age là một biến số liên tục. Để vẽ biểu đồ tần số của biến số age, chúng ta chỉ đơn gi n lệnh hist(age). Nh đã đề cập trên, chúng ta có thể c i tiến đồ thị này bằng cách cho thêm tựa đề chính (main) và tựa đề của trục hoành (xlab) và trục tung (ylab): > hist(age) > hist(age, main="Frequency distribution by age group", xlab="Age group", ylab="No of patients") Histogram of age 8 0 0 2 2 4 6 No of patients 6 4 Frequency 8 10 10 12 12 Frequency distribution by age group 40 50 60 70 80 age 40 50 60 70 80 Age group Biểu đồ 11a. Trục tung là số bệnh nhân (đối Biểu đồ 11b. Thêm tên biểu đồ và tên của trục t ợng nghiên cứu) và trục hoành là độ tuổi. trung và trục hoành bằng xlab và ylab. Chẳng h n nh tuổi 40 đến 45 có 6 bệnh nhân, từ 70 đến 80 tuổi có 4 bệnh nhân. Chúng ta cũng có thể biến đổi biểu đồ thành một đồ thị phân phối xác suất bằng hàm plot(density) nh sau (kết qu trong Biểu đồ 12a): > plot(density(age),add=TRUE) 48 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn density.default(x = age) Density 0.00 0.00 0.01 0.02 0.02 0.01 Density 0.03 0.03 0.04 0.04 Histogram of age 30 40 50 60 70 80 90 40 N = 50 Bandwidth = 3.806 50 60 70 80 age Biểu đồ 12a. Xác suất phân phối mật độ cho Biểu đồ 12b. Xác suất phân phối mật độ cho biến age (độ tuổi). biến age (độ tuổi) với nhiều interquartile. Chúng ta có thể vẽ hai đồ thị chồng lên bằng cách dùng hàm interquartile nh sau (kết qu xem Biểu đồ 12b): 8.6 Biểu đồ hộp (boxplot) Để vẽ biểu đồ hộp của biến số tc, chúng ta chỉ đơn gi n lệnh: > boxplot(tc, main="Box plot of total cholesterol", ylab="mg/L") 3 4 5 mg/L 6 7 8 Box plot of total cholesterol Biểu đồ 13. Trong biểu đồ này, chúng ta thấy median (trung vị) kho ng 5.6 mg/L, 25% total cholesterol thấp hơn 4.1, và 75% thấp hơn 6.2. Total cholesterol thấp nhất 49 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn là khoang 3, và cao nhất là trên 8 mg/L. Trong biểu đồ sau đây, chúng ta so sánh tc giữa hai nhóm nam và nữ: > boxplot(tc ~ sex, main=”Box plot of total cholestrol by sex”, ylab="mg/L") Kết qu trình bày trong Biểu đồ 14a. Chúng ta có thể biến đổ giao diện của đồ thị bằng cách dùng thông số horizontal=TRUE và thay đổi màu bằng thông số col nh sau (Biểu đồ 14b): > boxplot(tc~sex, horizontal=TRUE, main="Box plot of total cholesterol", ylab="mg/L", col = "pink") Box plot of total cholesterol 3 4 Nam 5 mg/L mg/L 6 7 Nu 8 Box plot of total cholesterol by sex 3 Nam 4 5 6 7 8 Nu Biểu đồ 14a. Trong biểu đồ này, chúng ta Biểu đồ 14b. Total cholesterol cho từng thấy trung vị của total cholesterol nữ giới giới tính, với màu sắc và hình hộp nằm thấp hơn nam giới, nh ng độ dao động giữa ngang. hai nhóm không khác nhau bao nhiêu. 8.7 Phân tích biểu đồ cho hai biến liên tục 8.7.1 Biểu đồ tán xạ (scatter plot) Để tìm hiểu mối liên hệ giữa hai biến, chúng ta dùng biểu đồ tán x . Để vẽ biểu đồ tán x về mối liên hệ giữa biến số tc và hdl, chúng ta sử dụng hàm plot. Thông số thứ nhất của hàm plot là trục hoành (x-axis) và thông số thứ 2 là trục tung. Để tìm hiểu mối liên hệ giữa tc và hdl chúng ta đơn gi n lệnh: > plot(tc, hdl) 50 Phân tích số liệu và biểu đồ bằng R 2 4 hdl 6 8 Nguyễn Văn Tuấn 3 4 5 6 7 8 tc Biểu đồ 15. Mối liên hệ giữa tc và hdl. Trong biểu đồ này, chúng ta vẽ biến số hdl trên trục tung và tc trên trục hoành. Chúng ta muốn phân biệt giới tính (nam và nữ) trong biểu đồ trên. Để vẽ biểu đồ đó, chúng ta ph i dùng đến hàm ifelse. Trong lệnh sau đây, nếu sex==”Nam” thì vẽ kí tự số 16 (ô tròn), nếu không nam thì vẽ kí tự số 22 (tức ô vuông): > plot(hdl, tc, pch=ifelse(sex=="Nam", 16, 22)) Kết qu là Biểu đồ 16a. Chúng ta cũng có thể thay kí t thành “M” (nam) và “F” nữ(xem Biểu đồ 16b): > plot(hdl, tc, pch=ifelse(sex=="Nam", “M”, “F”)) 51 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn M F 8 F 8 M F 7 F M 6 M F M M F tc hdl 6 M F F F M F M F M M F F 5 4 M M F F M F F F F 2 4 F M F M M F F 3 F 3 4 5 6 7 8 F 2 M M F 4 tc 6 8 hdl Biểu đồ 16a. Mối liên hệ giữa tc và hdl theo Biểu đồ 16a. Mối liên hệ giữa tc và hdl theo từng giới tính đ ợc thể hiện bằng hai kí hiệu từng giới tính đ ợc thể hiện bằng hai kí tự. dấu. Chúng ta cũng có thể vẽ một đ ng biểu diễn hồi qui tuyến tính (regression line) qua các điểm trên bằng cách tiếp tục ra các lệnh sau đây: > plot(hdl ~ tc, pch=16, main="Total cholesterol and HDL cholesterol", xlab="Total cholesterol", ylab="HDL cholesterol", bty=”l”) > reg <- lm(hdl ~ tc) > abline(reg) Kết qu là Biểu đồ 17a d ới đây. Chúng ta cũng có thể dùng hàm trơn (smooth function) để biểu diễn mối liên hệ giữa hai biến số. Đồ thị sau đây sử dụng lowess (một hàm thông th ng nhất) trong việc “làm trơn” số liệu tc và hdl (Biểu đồ 17b). > plot(hdl ~ tc, pch=16, main="Total cholesterol and HDL cholesterol with LOEWSS smooth function", xlab="Total cholesterol", ylab="HDL cholesterol", bty=”l”) > lines(lowess(hdl, tc, f=2/3, iter=3), col="red") 52 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn T otal cholesterol and HDL cholesterol 6 2 4 HDL cholesterol 4 2 HDL cholesterol 6 8 8 T otal cholesterol and HDL cholesterol with LOEWSS smooth function 3 4 5 6 7 8 3 Total cholesterol 4 5 6 7 8 Total cholesterol Biểu đồ 17a. Trong lệnh trên, reg<- Biểu đồ 17b. Thay vì dùng abline, chúng ta lm(hdl~tc) có nghĩa là tìm ph ơng trình dùng hàm lowess để thể hiện mối liên hệ giữa liên hệ giữa hdl và tc bằng linear model tc và hdl. (lm) và đ8 t kết qu vào đối t ợng reg. Lệnh thứ hai abline(reg) yêu cầu R vẽ đ ng thẳng từ ph ơng trình trong reg B n đọc có thể thí nghiệm với nhiều thông số f=1/2, f=2/5, hay thậm chí f=1/10 sẽ thấy đồ thị biến đổi một cách “thú vị”. 8.8 Phân tích Biểu đồ cho nhiều biến: pairs Chúng ta có thể tìm hiểu mối liên hệ giữa các biến số nh age, bmi, hdl, ldl và tc bằng cách dùng lệnh pairs. Nh ng tr ớc hết, chúng ta ph i đ a các biến số này vào một data.frame chỉ gồm những biến số có thể vẽ đ ợc, và sau đó sử dụng hàm pairs trong R. > lipid <- data.frame(age,bmi,hdl,ldl,tc) > pairs(lipid, pch=16) Kết qu sẽ là: 53 Phân tích số liệu và biểu đồ bằng R 20 22 24 1 2 3 4 5 6 70 80 18 Nguyễn Văn Tuấn 22 24 50 60 age 6 8 18 20 bmi 4 5 6 2 4 hdl 6 7 8 1 2 3 ldl 3 4 5 tc 50 60 70 80 2 4 6 8 3 4 5 6 7 8 8.9 Biểu đồ với sai số chuẩn (standard error) Trong biểu đồ sau đây, chúng ta có 5 nhóm (biến số x đ ợc mô phỏng chứ không ph i số liệu thật), và mỗi nhóm có giá trị trung bình mean, và độ tin cậy 95% (lcl và ucl). Thông th ng lcl=mean-1.96*SE và ucl = mean+1.96*SE (SE là sai số chuẩn). Chúng ta muốn vẽ biểu đồ cho 5 nhóm với sai số chuẩn đó. Các lệnh và hàm sau đây sẽ cần thiết: > > > > > > group <- c(1,2,3,4,5) mean <- c(1.1, 2.3, 3.0, 3.9, 5.1) lcl <- c(0.9, 1.8, 2.7, 3.8, 5.0) ucl <- c(1.3, 2.4, 3.5, 4.1, 5.3) plot(group, mean, ylim=range(c(lcl, ucl))) arrows(group, ucl, group, lcl, length=0.5, angle=90, code=3) 54 Phân tích số liệu và biểu đồ bằng R 3 1 2 mean 4 5 Nguyễn Văn Tuấn 1 2 3 4 5 group 9. Phân tích thống kê mô tả 9.1 Thống kê mô tả (descriptive statistics, summary) Để minh họa cho việc áp dụng R vào thống kê mô t , tôi sẽ sử dụng một dữ liệu nghiên cứu có tên là igfdata. Trong nghiên cứu này, ngoài các chỉ số liên quan đến giới tính, độ tuổi, trọng l ợng và chiều cao, chúng tôi đo l ng các hormone liên quan đến tình tr ng tăng tr ng nh igfi, igfbp3, als, và các markers liên quan đến sự chuyển hóa của x ơng pinp, ictp và pinp. Có 100 đối t ợng nghiên cứu. Dữ liệu này đ ợc chứa trong directory c:\works\stats. Tr ớc hết, chúng ta cần ph i nhập dữ liệu vào R với những lệnh sau đây (các câu chữ theo sau dấu # là những chú thích để b n đọc theo dõi): > options(width=100) # chuyển directory > setwd("c:/works/stats") # đọc dữ liệu vào R > igfdata <- read.table("igf.txt", header=TRUE, na.strings=".") > attach(igfdata) # xem xét các cột số trong dữ liệu > names(igfdata) [1] "id" "sex" "age" [7] "igfi" "igfbp3" "als" "weight" "pinp" "height" "ictp" "ethnicity" "p3np" > igfdata id sex age weight height ethnicity 55 igfi igfbp3 als pinp ictp p3np Phân tích số liệu và biểu đồ bằng R 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 ... ... 97 97 98 98 99 99 100 100 Nguyễn Văn Tuấn Female Male Female Female Female Female Female Female Female Female 15 16 15 15 16 25 19 18 15 24 42 44 43 42 47 45 45 43 41 45 162 Asian 189.000 4.00000 323.667 353.970 160 Caucasian 160.000 3.75000 333.750 375.885 157 Asian 146.833 3.43333 248.333 199.507 155 Asian 185.500 3.40000 251.000 483.607 167 Asian 192.333 4.23333 322.000 105.430 160 Asian 110.000 3.50000 284.667 76.487 161 Asian 157.000 3.20000 274.000 75.880 153 Asian 146.000 3.40000 303.000 86.360 149 Asian 197.667 3.56667 308.500 254.803 157 African 148.000 3.40000 273.000 44.720 11.2867 8.3367 10.4300 6.7450 8.3633 12.5000 13.3300 14.2767 7.9233 4.5033 4.9833 4.9367 6.3500 5.3200 7.3700 4.6700 11.8700 6.8200 3.7400 6.1600 Female Male Female Male 17 18 18 15 54 55 48 54 168 Caucasian 204.667 4.96667 441.333 64.130 5.1600 169 Asian 178.667 3.86667 273.000 185.913 7.5267 151 Asian 237.000 3.46667 324.333 105.127 5.9867 168 Asian 130.000 2.70000 259.333 325.840 10.2767 Trên đây chỉ là một phần số liệu trong số 100 đối t ợng. Cho một biến số x1 , x2 , x3 ,..., xn chúng ta có thể tính toán một số chỉ số thống kê mô t nh sau: Lí thuyết Số trung bình: x = Ph ơng sai: s 2 = 1 ∑ xi . n i =1 Hàm R mean(x) n 1 n 2 ∑ ( xi − x ) n − 1 i =1 var(x) Độ lệch chuẩn: s = s 2 Sai số chuẩn (standard error): SE = sd(x) Không có s n min(x) max(x) range(x) Trị số thấp nhất Trị số cao nhất Toàn cự (range) Ví dụ 9: Để tìm giá trị trung bình của độ tuổi, chúng ta chỉ đơn gi n lệnh: > mean(age) [1] 19.17 Hay ph ơng sai và độc lệch chuẩn của tuổi: > var(age) [1] 15.33444 > sd(age) [1] 3.915922 56 4.4367 8.8333 5.6600 6.5933 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Tuy nhiên, R có lệnh summary có thể cho chúng ta tất c thông tin thống kê về một biến số: > summary(age) Min. 1st Qu. 13.00 16.00 Median 19.00 Mean 3rd Qu. 19.17 21.25 Max. 34.00 Nói chung, kết qu này đơn gi n và các viết tắt cũng có thể dễ hiểu. Chú ý, trong kết qu trên, có hai chỉ số “1st Qu” và “3rd Qu” có nghĩa là first quartile (t ơng đ ơng với vị trí 25%) và third quartile (t ơng đ ơng với vị trí 75%) của một biến số. First quartile = 16 có nghĩa là 25% đối t ợng nghiên cứu có độ tuổi bằng hoặc nhỏ hơn 16 tuổi. T ơng tự, Third quartile = 34 có nghĩa là 75% đối t ợng có độ tuổi bằng hoặc thấp hơn 34 tuổi. Tất nhiên số trung vị (median) 19 cũng có nghĩa là 50% đối t ợng có độ tuổi 19 tr xuống (hay 19 tuổi tr lên). R không có hàm tính sai số chuẩn, và trong hàm summary, R cũng không cung cấp độ lệch chuẩn. Để có các số này, chúng ta có thể tự viết một hàm đơn gi n (hãy gọi là desc) nh sau: desc <- function(x) { av <- mean(x) sd <- sd(x) se <- sd/sqrt(length(x)) c(MEAN=av, SD=sd, SE=se) } Và có thể gọi hàm này để tính bất cứ biến nào chúng ta muốn, nh tính biến als sau đây: > desc(als) MEAN SD 301.841120 58.987189 SE 5.898719 Để có một “quang c nh” chung về dữ liệu igfdata chúng ta chỉ đơn gi n lệnh summary nh sau: > summary(igfdata) id sex Min. : 1.00 Female:69 1st Qu.: 25.75 Male :31 Median : 50.50 Mean : 50.50 3rd Qu.: 75.25 Max. :100.00 igfi igfbp3 age Min. :13.00 1st Qu.:16.00 Median :19.00 Mean :19.17 3rd Qu.:21.25 Max. :34.00 als 57 weight Min. :41.00 1st Qu.:47.00 Median :50.00 Mean :49.91 3rd Qu.:53.00 Max. :60.00 pinp height Min. :149.0 1st Qu.:157.0 Median :162.0 Mean :163.1 3rd Qu.:168.0 Max. :196.0 ictp ethnicity African : 8 Asian :60 Caucasian:30 Others : 2 Phân tích số liệu và biểu đồ bằng R Min. : 85.71 1st Qu.:137.17 Median :161.50 Mean :165.59 3rd Qu.:186.46 Max. :427.00 Min. :2.000 1st Qu.:3.292 Median :3.550 Mean :3.617 3rd Qu.:3.875 Max. :5.233 Nguyễn Văn Tuấn Min. :192.7 1st Qu.:256.8 Median :292.5 Mean :301.8 3rd Qu.:331.2 Max. :471.7 Min. : 26.74 1st Qu.: 68.10 Median :103.26 Mean :167.17 3rd Qu.:196.45 Max. :742.68 Min. : 2.697 1st Qu.: 4.878 Median : 6.338 Mean : 7.420 3rd Qu.: 8.423 Max. :21.237 p3np Min. : 2.343 1st Qu.: 4.433 Median : 5.445 Mean : 6.341 3rd Qu.: 7.150 Max. :16.303 R tính toán tất c các biến số nào có thể tính toán đ ợc! Thành ra, ngay c cột id (tức mã số của đối t ợng nghiên cứu) R cũng tính luôn! (và chúng ta biết kết qu của cột id chẳng có ý nghĩa thống kê gì). Đối với các biến số mang tính phân lo i nh sex và ethnicity (sắc tộc) thì R chỉ báo cáo tần số cho mỗi nhóm. Kết qu trên cho tất c đối t ợng nghiên cứu. Nếu chúng ta muốn kết qu cho từng nhóm nam và nữ riêng biệt, hàm by trong R rất hữu dụng. Trong lệnh sau đây, chúng ta yêu cầu R tóm l ợc dữ liệu igfdata theo sex. > by(igfdata, sex, summary) sex: Female id Min. : 1.0 1st Qu.:21.0 Median :47.0 Mean :48.2 3rd Qu.:75.0 Max. :99.0 ethnicity African : 4 Asian :43 Caucasian:22 Others : 0 sex Female:69 Male : 0 age weight height Min. :13.00 Min. :41.00 Min. :149.0 1st Qu.:17.00 1st Qu.:47.00 1st Qu.:156.0 Median :19.00 Median :50.00 Median :162.0 Mean :19.59 Mean :49.35 Mean :161.9 3rd Qu.:22.00 3rd Qu.:52.00 3rd Qu.:166.0 Max. :34.00 Max. :60.00 Max. :196.0 igfi igfbp3 als Min. : 85.71 Min. :2.767 Min. :204.3 1st Qu.:136.67 1st Qu.:3.333 1st Qu.:263.8 Median :163.33 Median :3.567 Median :302.7 Mean :167.97 Mean :3.695 Mean :311.5 3rd Qu.:186.17 3rd Qu.:3.933 3rd Qu.:361.7 Max. :427.00 Max. :5.233 Max. :471.7 pinp ictp p3np Min. : 26.74 Min. : 2.697 Min. : 2.343 1st Qu.: 62.75 1st Qu.: 4.717 1st Qu.: 4.337 Median : 78.50 Median : 5.537 Median : 5.143 Mean :108.74 Mean : 6.183 Mean : 5.643 3rd Qu.:115.26 3rd Qu.: 7.320 3rd Qu.: 6.143 Max. :502.05 Max. :13.633 Max. :14.420 -----------------------------------------------------------sex: Male id sex age weight height Min. : 2.00 Female: 0 Min. :14.00 Min. :44.00 Min. :155.0 1st Qu.: 34.50 Male :31 1st Qu.:15.00 1st Qu.:48.50 1st Qu.:161.5 Median : 56.00 Median :17.00 Median :51.00 Median :164.0 Mean : 55.61 Mean :18.23 Mean :51.16 Mean :165.6 3rd Qu.: 75.00 3rd Qu.:20.00 3rd Qu.:53.50 3rd Qu.:169.0 Max. :100.00 Max. :27.00 Max. :59.00 Max. :191.0 ethnicity igfi igfbp3 als 58 Phân tích số liệu và biểu đồ bằng R African : 4 Asian :17 Caucasian: 8 Others : 2 pinp Min. : 56.28 1st Qu.:135.07 Median :245.92 Mean :297.21 3rd Qu.:450.38 Max. :742.68 Min. : 94.67 1st Qu.:138.67 Median :160.00 Mean :160.29 3rd Qu.:183.00 Max. :274.00 ictp Min. : 3.650 1st Qu.: 6.900 Median : 9.513 Mean :10.173 3rd Qu.:13.517 Max. :21.237 Nguyễn Văn Tuấn Min. :2.000 Min. :192.7 1st Qu.:3.183 1st Qu.:249.8 Median :3.500 Median :276.0 Mean :3.443 Mean :280.2 3rd Qu.:3.775 3rd Qu.:311.3 Max. :4.500 Max. :388.7 p3np Min. : 3.390 1st Qu.: 5.375 Median : 7.140 Mean : 7.895 3rd Qu.:10.010 Max. :16.303 Để xem qua phân phối của các hormones và chỉ số sinh hóa cùng một lúc, chúng ta có thể vẽ đồ thị cho tất c 6 biến số. Tr ớc hết, chia màn nh thành 6 cửa sổ (với 2 dòng và 3 cột); sau đó lần l ợt vẽ: > > > > > > > op <- par(mfrow=c(2,3)) hist(igfi) hist(igfbp3) hist(als) hist(pinp) hist(ictp) hist(p3np) 59 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Histogram of igfbp3 Histogram of als 200 300 400 0 0 100 20 Frequency 10 20 Frequency 10 20 0 10 Frequency 30 30 30 40 40 Histogram of igfi 2.0 3.0 4.0 5.0 150 250 350 450 igfbp3 als Histogram of pinp Histogram of ictp Histogram of p3np 40 30 20 Frequency 30 0 200 400 pinp 600 800 0 0 0 10 10 10 20 Frequency 30 20 Frequency 40 50 igf i 5 10 15 20 ictp 5 10 15 p3np 9.2 Thống kê mô tả theo từng nhóm Nếu chúng ta muốn tính trung bình của một biến số nh igfi cho mỗi nhóm nam và nữ giới, hàm tapply trong R có thể dùng cho việc này: > tapply(igfi, list(sex), mean) Female Male 167.9741 160.2903 Trong lệnh trên, igfi là biến số chúng ta cần tính, biến số phân nhóm là sex, và chỉ số thống kê chúng ta muốn là trung bình (mean). Qua kết qu trên, chúng ta thấy số trung bình của igfi cho nữ giới (167.97) cao hơn nam giới (160.29). Nh ng nếu chúng ta muốn tính cho từng giới tính và sắc tộc, chúng ta chỉ cần thêm một biến số trong hàm list: > tapply(igfi, list(ethnicity, sex), mean) Female Male African 145.1252 120.9168 60 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Asian 165.6589 160.4999 Caucasian 176.6536 169.4790 Others NA 200.5000 Trong kết qu trên, NA có nghĩa là “not available”, tức không có số liệu cho phụ nữ trong các sắc tộc “others”. 9.3 Kiểm định t (t.test) Kiểm định t dựa vào gi thiết phân phối chuẩn. Có hai lo i kiểm định t: kiểm định t cho một mẫu (one-sample t-test), và kiểm định t cho hai mẫu (two-sample t-test). Kiểm định t một mẫu nằm tr l i câu hỏi dữ liệu từ một mẫu có ph i thật sự bằng một thông số nào đó hay không. Còn kiểm định t hai mẫu thì nhằm tr l i câu hỏi hai mẫu có cùng một luật phân phối, hay cụ thể hơn là hai mẫu có thật sự có cùng trị số trung bình hay không. Tôi sẽ lần l ợt minh họa hai kiểm định này qua số liệu igfdata trên. 9.3.1 Kiểm định t một mẫu Ví dụ 10. Qua phân tích trên, chúng ta thấy tuổi trung bình của 100 đối t ợng trong nghiên cứu này là 19.17 tuổi. Chẳng h n nh trong quần thể này, tr ớc đây chúng ta biết rằng tuổi trung bình là 30 tuổi. Vấn đề đặt ra là có ph i mẫu mà chúng ta có đ ợc có đ i diện cho quần thể hay không. Nói cách khác, chúng ta muốn biết giá trị trung bình 19.17 có thật sự khác với giá trị trung bình 30 hay không. Để tr l i câu hỏi này, chúng ta sử dụng kiểm định t. Theo lí thuyết thống kê, kiểm định t đ ợc định nghĩa bằng công thức sau đây: t= x −µ s/ n Trong đó, x là giá trị trung bình của mẫu, µ là trung bình theo gi thiết (trong tr ng hợp này, 30), s là độ lệch chuẩn, và n là số l ợng mẫu (100). Nếu giá trị t cao hơn giá trị lí thuyết theo phân phối t một tiêu chuẩn có ý nghĩa nh 5% chẳng h n thì chúng ta có lí do để phát biểu khác biệt có ý nghĩa thống kê. Giá trị này cho mẫu 100 có thể tính toán bằng hàm qt của R nh sau: > qt(0.95, 100) [1] 1.660234 Nh ng có một cách tính toán nhanh gọn hơn để tr l i câu hỏi trên, bằng cách dùng hàm t.test nh sau: > t.test(age, mu=30) One Sample t-test 61 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn data: age t = -27.6563, df = 99, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 30 95 percent confidence interval: 18.39300 19.94700 sample estimates: mean of x 19.17 Trong lệnh trên age là biến số chúng ta cần kiểm định, và mu=30 là giá trị gi thiết. R trình bày trị số t = -27.66, với 99 bậc tự do, và trị số p < 2.2e-16 (tức rất thấp). R cũng cho biết độ tin cậy 95% của age là từ 18.4 tuổi đến 19.9 tuổi (30 tuổi nằm quá ngoài kho ng tin cậy này). Nói cách khác, chúng ta có lí do để phát biểu rằng độ tuổi trung bình trong mẫu này thật sự thấp hơn độ tuổi trung bình của quần thể. 9.3.2 Kiểm định t hai mẫu Ví dụ 11. Qua phân tích mô t trên (phầm summary) chúng ta thấy phụ nữ có độ hormone igfi cao hơn nam giới (167.97 và 160.29). Câu hỏi đặt ra là có ph i thật sự đó là một khác biệt có hệ thống hay do các yếu tố ngẫu nhiên gây nên. Tr l i câu hỏi này, chúng ta cần xem xét mức độ khác biệt trung bình giữa hai nhóm và độ lệch chuẩn của độ khác biệt. x2 − x1 SED Trong đó x1 và x2 là số trung bình của hai nhóm nam và nữ, và SED là độ lệch chuẩn của ( x1 - x2 ) . Thực ra, SED có thể ớc tính bằng công thức: t= SED = SE12 + SE22 Trong đó SE1 và SE2 là sai số chuẩn (standard error) của hai nhóm nam và nữ. Theo lí thuyết xác suất, t tuân theo luật phân phối t với bậc tự do n1 + n2 − 2 , trong đó n1 và n2 là số mẫu của hai nhóm. Chúng ta có thể dùng R để tr l i câu hỏi trên bằng hàm t.test nh sau: > t.test(igfi~ sex) Welch Two Sample t-test data: igfi by sex t = 0.8412, df = 88.329, p-value = 0.4025 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -10.46855 25.83627 sample estimates: mean in group Female mean in group Male 167.9741 160.2903 62 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn R trình bày các giá trị quan trọng tr ớc hết: t = 0.8412, df = 88.329, p-value = 0.4025 df là bậc tự do. Trị số p = 0.4025 cho thấy mức độ khác biệt giữa hai nhóm nam và nữ không có ý nghĩa thống kê (vì cao hơn 0.05 hay 5%). 95 percent confidence interval: -10.46855 25.83627 là kho ng tin cậy 95% về độ khác biệt giữa hai nhóm. Kết qu tính toán trên cho biết độ igf nữ giới có thể thấp hơn nam giới 10.5 ng/L hoặc cao hơn nam giới kho ng 25.8 ng/L. Vì độ khác biệt quá lớn và đó là thêm bằng chứng cho thấy không có khác biệt có ý nghĩa thống kê giữa hai nhóm. Kiểm định trên dựa vào gi thiết hai nhóm nam và nữ có khác ph ơng sai. Nếu chúng ta có lí do đề cho rằng hai nhóm có cùng ph ơng sai, chúng ta chỉ thay đổi một thông số trong hàm t với var.equal=TRUE nh sau: > t.test(igfi~ sex, var.equal=TRUE) Two Sample t-test data: igfi by sex t = 0.7071, df = 98, p-value = 0.4812 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -13.88137 29.24909 sample estimates: mean in group Female mean in group Male 167.9741 160.2903 Về mặc số, kết qu phân tích trên có khác chút ít so với kết qu phân tích dựa vào gi định hai ph ơng sai khác nhau, nh ng trị số p cũng đi đến một kết luận rằng độ khác biệt giữa hai nhóm không có ý nghĩa thống kê. 9.4 Kiểm định Wilcoxon cho hai mẫu (wilcox.test) Kiểm định t dựa vào gi thiết là phân phối của một biến ph i tuân theo luật phân phối chuẩn. Nếu gi định này không đúng, kết qu của kiểm định t có thể không hợp lí (valid). Để kiểm định phân phối của igfi, chúng ta có thể dùng hàm shapiro.test nh sau: > shapiro.test(igfi) Shapiro-Wilk normality test 63 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn data: igfi W = 0.8528, p-value = 1.504e-08 Trị số p nhỏ hơn 0.05 rất nhiều, cho nên chúng ta có thể nói rằng phân phối của igfi không tuân theo luật phân phối chuẩn. Trong tr ng hợp này, việc so sánh giữa hai nhóm có thể dựa vào ph ơng pháp phi tham số (non-parametric) có tên là kiểm định Wilcoxon, vì kiểm định này (không nh kiểm định t) không tùy thuộc vào gi định phân phối chuẩn. > wilcox.test(igfi ~ sex) Wilcoxon rank sum test with continuity correction data: igfi by sex W = 1125, p-value = 0.6819 alternative hypothesis: true mu is not equal to 0 Trị số p = 0.682 cho thấy qu thật độ khác biệt về igfi giữa hai nhóm nam và nữ không có ý nghĩa thống kê. Kết luận này cũng không khác với kết qu phân tích bằng kiểm định t. 9.5 Kiểm định t cho các biến số theo cặp (paired t-test, t.test) Kiểm định t vừa trình bày trên là cho các nghiên cứu gồm hai nhóm độc lập nhau (nh giữa hai nhóm nam và nữ), nh ng không thể ứng dụng cho các nghiên cứu mà một nhóm đối t ợng đ ợc theo dõi theo th i gian. Tôi t m gọi các nghiên cứu này là nghiên cứu theo cặp. Trong các nghiên cứu này, chúng ta cần sử dụng một kiểm định t có tên là paired t-test. Ví dụ 12. Một nhóm bệnh nhân gồm 10 ng i đ ợc điều trị bằng một thuốc nhằm gi m huyết áp. Huyết áp của bệnh nhân đ ợc đo lúc kh i đầu nghiên cứu (lúc ch a điều trị), và sau khi điều khị. Số liệu huyết áp của 10 bệnh nhân nh sau: Tr ớc khi điều trị (x0) Sau khi điều trị (x1) 180, 140, 160, 160, 220, 185, 145, 160, 160, 170 170, 145, 145, 125, 205, 185, 150, 150, 145, 155 Câu hỏi đặt ra là độ biến chuyển huyết áp trên có đủ để kết luận rằng thuốc điều trị có hiệu qu gi m áp huyết. Để tr l i câu hỏi này, chúng ta dùng kiểm định t cho từng cặp nh sau: > > > > # nhập dữ kiện before <- c(180, 140, 160, 160, 220, 185, 145, 160, 160, 170) after <- c(170, 145, 145, 125, 205, 185, 150, 150, 145, 155) bp <- data.frame(before, after) > # kiểm định t > t.test(before, after, paired=TRUE) 64 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Paired t-test data: before and after t = 2.7924, df = 9, p-value = 0.02097 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 1.993901 19.006099 sample estimates: mean of the differences 10.5 Kết qu trên cho thấy sau khi điều trị áp suất máu gi m 10.5 mmHg, và kho ng tin cậy 95% là từ 2.0 mmHg đến 19 mmHg, với trị số p = 0.0209. Nh vậy, chúng ta có bằng chứng để phát biểu rằng mức độ gi m huyết áp có ý nghĩa thống kê. Chú ý nếu chúng ta phân tích sai bằng kiểm định thống kê cho hai nhóm độc lập d ới đây thì trị số p = 0.32 cho biết mức độ gi m áp suất không có ý nghĩa thống kê! > t.test(before, after) Welch Two Sample t-test data: before and after t = 1.0208, df = 17.998, p-value = 0.3209 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -11.11065 32.11065 sample estimates: mean of x mean of y 168.0 157.5 9.6 Kiểm định Wilcoxon cho các biến số theo cặp (wilcox.test) Thay vì dùng kiểm định t cho từng cặp, chúng ta cũng có thể sử dụng hàm wilcox.test cho cùng mục đích: > wilcox.test(before, after, paired=TRUE) Wilcoxon signed rank test with continuity correction data: before and after V = 42, p-value = 0.02291 alternative hypothesis: true mu is not equal to 0 Kết qu trên một lần nữa khẳng định rằng độ gi m áp suất máu có ý nghĩa thống kê với trị số (p=0.023) chẳng khác mấy so với kiểm định t cho từng cặp. 65 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 9.7 Tần số (frequency) Hàm table trong R có chức năng cho chúng ta biết về tần số của một biến số mang tính phân lo i nh sex và ethnicity. > table(sex) sex Female Male 69 31 > table(ethnicity) ethnicity African Asian Caucasian 8 60 30 Others 2 Một b ng thống kê 2 chiều: > table(sex, ethnicity) ethnicity sex African Asian Caucasian Others Female 4 43 22 0 Male 4 17 8 2 Chú ý trong các b ng thống kê trên, hàm table không cung cấp cho chúng ta số phần trăm. Để tính số phần trăm, chúng ta cần đến hàm prop.table và cách sử dụng có thể minh ho nh sau: # tạo ra một object tên là freq để chứa kết quả tần số > freq <- table(sex, ethnicity) # kiểm tra kết quả > freq ethnicity sex African Asian Caucasian Others Female 4 43 22 0 Male 4 17 8 2 # dùng hàm margin.table để xem kết quả > margin.table(freq, 1) sex Female Male 69 31 > margin.table(freq, 2) ethnicity African Asian Caucasian Others 66 Phân tích số liệu và biểu đồ bằng R 8 60 Nguyễn Văn Tuấn 30 2 # tính phần trăm bằng hàm prop.table > prop.table(freq, 1) ethnicity sex African Asian Caucasian Others Female 0.05797101 0.62318841 0.31884058 0.00000000 Male 0.12903226 0.54838710 0.25806452 0.06451613 Trong b ng thống kê trên, prop.table tính tỉ lệ sắc tộc cho từng giới tính. Chẳng h n nh nữ giới (female), 5.8% là ng i Phi châu, 62.3% là ng i Á châu, 31.8% là ng i Tây ph ơng da trắng . Tổng cộng là 100%. T ơng tự, ớ nam giới tỉ lệ ng i Phi châu là 12.9%, Á châu là 54.8%, v.v… # tính phần trăm bằng hàm prop.table > prop.table(freq, 2) ethnicity sex African Asian Caucasian Others Female 0.5000000 0.7166667 0.7333333 0.0000000 Male 0.5000000 0.2833333 0.2666667 1.0000000 Trong b ng thống kê trên, prop.table tính tỉ lệ giới tính cho từng sắc tộc. Chẳng h n nh trong nhóm ng i Á châu, 71.7% là nữ và 28.3% là nam. # tính phần trăm cho toàn bộ bảng > freq/sum(freq) ethnicity sex African Asian Caucasian Others Female 0.04 0.43 0.22 0.00 Male 0.04 0.17 0.08 0.02 9.8 Kiểm định tỉ lệ (proportion test, prop.test, binom.test) Kiểm định một tỉ lệ th ng dựa vào gi định phân phối nhị phân (binomial distribution). Với một số mẫu n và tỉ lệ p, và nếu n lớn (tức hơn 50 chẳng h n), thì phân phối nhị phân có thể t ơng đ ơng với phân phối chuẩn với số trung bình np và ph ơng sai np(1 – p). Gọi x là số biến cố mà chúng ta quan tâm, kiểm định gi thiết p = π có thể sử dụng thống kê sau đây: z= x − nπ nπ (1 − π ) đây, z tuân theo luật phân phối chuẩn với trung bình 0 và ph ơng sai 1. Cũng có thể nói z2 tuân theo luật phân phối Chi bình ph ơng với bậc tự do bằng 1. 67 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ví dụ 13. Trong nghiên cứu trên, chúng ta thấy có 69 nữ và 31 nam. Nh vậy tỉ lệ nữ là 0.69 (hay 69%). Để kiểm định xem tỉ lệ này có thật sự khác với tỉ lệ 0.5 hay không, chúng ta có thể sử dụng hàm prop.test(x, n, π) nh sau: > prop.test(69, 100, 0.50) 1-sample proportions test with continuity correction data: 69 out of 100, null probability 0.5 X-squared = 13.69, df = 1, p-value = 0.0002156 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.5885509 0.7766330 sample estimates: p 0.69 Trong kết qu trên, prop.test ớc tính tỉ lệ nữ giới là 0.69, và kho ng tin cậy 95% là 0.588 đến 0.776. Giá trị Chi bình ph ơng là 13.69, với trị số p = 0.00216. Nh vậy, nghiên cứu này có tỉ lệ nữ cao hơn 50%. Một cách tính chính xác hơn kiểm định tỉ lệ là kiểm định nhị phân bionom.test(x, n, π) nh sau: > binom.test(69, 100, 0.50) Exact binomial test data: 69 and 100 number of successes = 69, number of trials = 100, p-value = 0.0001831 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.5896854 0.7787112 sample estimates: probability of success 0.69 Nói chung, kết qu của kiểm định nhị phân không khác gì so với kiểm định Chi bình ph ơng, với trị số p = 0.00018, chúng ta càng có bằng chứng để kết luận rằng tỉ lệ nữ giới trong nghiên cứu này thật sự cao hơn 50%. 9.9 So sánh hai tỉ lệ (prop.test, binom.test) Ph ơng pháp so sánh hai tỉ lệ có thể khai triển trực tiếp từ lí thuyết kiểm định một tỉ lệ vừa trình bày trên. Cho hai mẫu với số đối t ợng n1 và n2, và số biến cố là x1 và x2. Do đó, chúng ta có thể ớc tính hai tỉ lệ p1 và p2. Lí thuyết xác suất cho phép chúng ta phát biểu rằng độ khác biệt giữa hai mẫu d = p1 – p2 tuân theo luật phân phối chuẩn với số trung bình 0 và ph ơng sai bằng: 68 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 1 1 Vd =  +  p (1 − p )  n1 n2  Trong đó: p= x1 + x2 n1 + n2 Thành ra, z = d/Vd tuân theo luật phân phối chuẩn với trung bình 0 và ph ơng sai 1. Nói cách khác, z2 tuân theo luật phân phối Chi bình ph ơng với bậc tự do bằng 1. Do đó, chúng ta cũng có thể sử dụng prop.test để kiểm định hai tỉ lệ. Ví dụ 14. Một nghiên cứu đ ợc tiến hành so sánh hiệu qu của thuốc chống gãy x ơng. Bệnh nhân đ ợc chia thành hai nhóm: nhóm A đ ợc điều trị gồm có 100 bệnh nhân, và nhóm B không đ ợc điều trị gồm 110 bệnh nhân. Sau th i gian 12 tháng theo dõi, nhóm A có 7 ng i bị gãy x ơng, và nhóm B có 20 ng i gãy x ơng. Vấn đề đặt ra là tỉ lệ gãy x ơng trong hai nhóm này bằng nhau (tức thuốc không có hiệu qu )? Để kiểm định xem hai tỉ lệ này có thật sự khác nhau, chúng ta có thể sử dụng hàm prop.test(x, n, π) nh sau: > fracture <- c(7, 20) > total <- c(100, 110) > prop.test(fracture, total) 2-sample test for equality of proportions with continuity correction data: fracture out of total X-squared = 4.8901, df = 1, p-value = 0.02701 alternative hypothesis: two.sided 95 percent confidence interval: -0.20908963 -0.01454673 sample estimates: prop 1 prop 2 0.0700000 0.1818182 Kết qu phân tích trên cho thấy tỉ lệ gãy x ơng trong nhóm 1 là 0.07 và nhóm 2 là 0.18. Phân tích trên còn cho thấy xác suất 95% rằng độ khác biệt giữa hai nhóm có thể 0.01 đến 0.20 (tức 1 đến 20%). Với trị số p = 0.027, chúng ta có thể nói rằng tỉ lệ gãy x ơng trong nhóm A qu thật thấp hơn nhóm B. 9.10 So sánh nhiều tỉ lệ (prop.test, chisq.test) Kiểm định prop.test còn có thể sử dụng để kiểm định nhiều tỉ lệ cùng một lúc. Trong nghiên cứu trên, chúng ta có 4 nhóm sắc tộc và tần số cho từng giới tính nh sau: > table(sex, ethnicity) 69 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn ethnicity sex African Asian Caucasian Others Female 4 43 22 0 Male 4 17 8 2 Chúng ta muốn biết tỉ lệ nữ giới giữa 4 nhóm sắc tộc có khác nhau hay không, và để tr l i câu hỏi này, chúng ta l i dùng prop.test nh sau: > female <- c( 4, 43, 22, 0) > total <- c(8, 60, 30, 2) > prop.test(female, total) 4-sample test for equality of proportions without continuity correction data: female out of total X-squared = 6.2646, df = 3, p-value = 0.09942 alternative hypothesis: two.sided sample estimates: prop 1 prop 2 prop 3 prop 4 0.5000000 0.7166667 0.7333333 0.0000000 Warning message: Chi-squared approximation may be incorrect in: prop.test(female, total) Tuy tỉ lệ nữ giới giữa các nhóm có vẻ khác nhau lớn (73% trong nhóm 3 (ng i da trắng) so với 50% trong nhóm 1 (Phi châu) và 71.7% trong nhóm Á châu, nh ng kiểm định Chi bình ph ơng cho biết trên ph ơng diện thống kê, các tỉ lệ này không khác nhau, vì trị số p = 0.099. 9.10.1 Kiểm định Chi bình ph ơng (Chi squared test, chisq.test) Thật ra, kiểm định Chi bình ph ơng còn có thể tính toán bằng hàm chisq.test nh sau: > chisq.test(sex, ethnicity) Pearson's Chi-squared test data: sex and ethnicity X-squared = 6.2646, df = 3, p-value = 0.09942 Warning message: Chi-squared approximation ethnicity) may be incorrect Kết qu này hoàn toàn giống với kết qu từ hàm prop.test. 70 in: chisq.test(sex, Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 9.10.2 Kiểm định Fisher (Fisher’s exact test, fisher.test) Trong kiểm định Chi bình ph ơng trên, chúng ta chú ý c nh báo: “Warning message: Chi-squared approximation may be incorrect in: prop.test(female, total)” Vì trong nhóm 4, không có nữ giới cho nên tỉ lệ là 0%. Hơn nữa, trong nhóm này chỉ có 2 đối t ợng. Vì số l ợng đối t ợng quá nhỏ, cho nên các ớc tính thống kê có thể không đáng tin cậy. Một ph ơng pháp khác có thể áp dụng cho các nghiên cứu với tần số thấp nh trên là kiểm định fisher (còn gọi là Fisher’s exact test). B n đọc có thể tham kh o lí thuyết đằng sau kiểm định fisher để hiểu rõ hơn về logic của ph ơng pháp này, nh ng đây, chúng ta chỉ quan tâm đến cách dùng R để tính toán kiểm định này. Chúng ta chỉ đơn gi n lệnh: > fisher.test(sex, ethnicity) Fisher's Exact Test for Count Data data: sex and ethnicity p-value = 0.1048 alternative hypothesis: two.sided Chú ý trị số p từ kiểm định Fisher là 0.1048, tức rất gần với trị số p của kiểm định Chi bình ph ơng. Cho nên, chúng ta có thêm bằng chứng để khẳng định rằng tỉ lệ nữ giới giữa các sắc tộc không khác nhau một cách đáng kể. 10. Phân tích hồi qui tuyến tính Ví dụ 15. Để minh họa cho vấn đề, chúng ta thử xem xét nghiên cứu sau đây, mà trong đó nhà nghiên cứu đo l ng độ cholestrol trong máu của 18 đối t ợng nam. Tỉ trọng cơ thể (body mass index) cũng đ ợc ớc tính cho mỗi đối t ợng bằng công thức tính BMI là lấy trọng l ợng (tính bằng kg) chia cho chiều cao bình ph ơng (m2). Kết qu đo l ng nh sau: Độ tuổi, tỉ trọng cơ thể và cholesterol Mã số ID (id) 1 2 3 4 5 6 7 8 Độ tuổi (age) 46 20 52 30 57 25 28 36 BMI (bmi) 25.4 20.6 26.2 22.6 25.4 23.1 22.7 24.9 Cholesterol (chol) 3.5 1.9 4.0 2.6 4.5 3.0 2.9 3.8 71 Phân tích số liệu và biểu đồ bằng R 9 10 11 12 13 14 15 16 17 18 Nguyễn Văn Tuấn 22 43 57 33 22 63 40 48 28 49 19.8 25.3 23.2 21.8 20.9 26.7 26.4 21.2 21.2 22.8 2.1 3.8 4.1 3.0 2.5 4.6 3.2 4.2 2.3 4.0 Nhìn sơ qua số liệu chúng ta thấy ng i có độ tuổi càng cao độ cholesterol cũng càng cao. Chúng ta thử nhập số liệu này vào R và vẽ một biểu đồ tán x nh sau: > age <- c(46,20,52,30,57,25,28,36,22,43,57,33,22,63,40,48,28,49) > bmi <-c(25.4,20.6,26.2,22.6,25.4,23.1,22.7,24.9,19.8,25.3,23.2, 21.8,20.9,26.7,26.4,21.2,21.2,22.8) > chol <- c(3.5,1.9,4.0,2.6,4.5,3.0,2.9,3.8,2.1,3.8,4.1,3.0, 2.5,4.6,3.2, 4.2,2.3,4.0) 2.0 2.5 3.0 chol 3.5 4.0 4.5 > data <- data.frame(age, bmi, chol) > plot(chol ~ age, pch=16) 20 30 40 50 60 age Biểu đồ 18. Liên hệ giữa độ tuổi và cholesterol. 72 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Biểu đồ 18 trên đây gợi ý cho thấy mối liên hệ giữa độ tuổi (age) và cholesterol là một đ ng thẳng (tuyến tính). Để “đo l ng” mối liên hệ này, chúng ta có thể sử dụng hệ số t ơng quan (coefficient of correlation). 10.1 Hệ số t ơng quan Hệ số t ơng quan (r) là một chỉ số thống kê đo l ng mối liên hệ t ơng quan giữa hai biến số, nh giữa độ tuổi (x) và cholesterol (y). Hệ số t ơng quan có giá trị từ -1 đến 1. Hệ số t ơng quan bằng 0 (hay gần 0) có nghĩa là hai biến số không có liên hệ gì với nhau; ng ợc l i nếu hệ số bằng -1 hay 1 có nghĩa là hai biến số có một mối liên hệ tuyệt đối. Nếu giá trị của hệ số t ơng quan là âm (r <0) có nghĩa là khi x tăng cao thì y gi m (và ng ợc l i, khi x gi m thì y tăng); nếu giá trị hệ số t ơng quan là d ơng (r > 0) có nghĩa là khi x tăng cao thì y cũng tăng, và khi x tăng cao thì y cũng gi m theo. Thực ra có nhiều hệ số t ơng quan trong thống kê, nh ng đây tôi sẽ trình bày 3 hệ số t ơng quan thông dụng nhất: hệ số t ơng quan Pearson r, Spearman ρ, và Kendall τ. 10.1.1 Hệ số t ơng quan Pearson Cho hai biến số x và y từ n mẫu, hệ số t ơng quan Pearson đ ợc ớc tính bằng công thức sau đây: r = ∑ ( xi − x )( yi − y ) n i =1 ∑ ( xi − x ) ∑ ( yi − y ) n i =1 2 n . Trong đó, nh định nghĩa phần trên, x 2 i =1 và y là giá trị trung bình của biến số x và y. Để ớc tính hệ số t ơng quan giữa độ tuổi age và cholesterol, chúng ta có thể sử dụng hàm cor(x,y) nh sau: > cor(age, chol) [1] 0.936726 Chúng ta có thể kiểm định gi thiết hệ số t ơng quan bằng 0 (tức hai biến x và y không có liên hệ). Ph ơng pháp kiểm định này th ng dựa vào phép biến đổi Fisher mà R đã có sẵn một hàm cor.test để tiến hành việc tính toán. > cor.test(age, chol) Pearson's product-moment correlation data: age and chol t = 10.7035, df = 16, p-value = 1.058e-08 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.8350463 0.9765306 sample estimates: cor 0.936726 73 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 10.1.2 Hệ số t ơng quan Spearman ρ Hệ số t ơng quan Pearson chỉ hợp lí nếu biến số x và y tuân theo luật phân phối chuẩn. Nếu x và y không tuân theo luật phân phối chuẩn, chúng ta ph i sử dụng một hệ số t ơng quan khác tên là Spearman, một ph ơng pháp phân tích phi tham số. Hệ số này đ ợc ớc tính bằng cách biến đổi hai biến số x và y thành thứ bậc (rank), và xem độ t ơng quan giữa hai dãy số bậc. Do đó, hệ số còn có tên tiếng Anh là Spearman’s Rank correlation. R ớc tính hệ số t ơng quan Spearman bằng hàm cor.test với thông số method=”spearman” nh sau: > cor.test(age, chol, method="spearman") Spearman's rank correlation rho data: age and chol S = 51.1584, p-value = 2.57e-09 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.947205 Warning message: Cannot compute exact p-values with ties in: cor.test.default(age, chol, method = "spearman") 10.1.3 Hệ số t ơng quan Kendall τ Hệ số t ơng quan Kendall (cũng là một ph ơng pháp phân tích phi tham số) đ ợc ớc tính bằng cách tìm các cặp số (x, y) “song hành" với nhau. Một cặp (x, y) song hành đây đ ợc định nghĩa là hiệu (độ khác biệt) trên trục hoành có cùng dấu hiệu (d ơng hay âm) với hiệu trên trục tung. Nếu hai biến số x và y không có liên hệ với nhau, thì số cặp song hành bằng hay t ơng đ ơng với số cặp không song hành. B i vì có nhiều cặp ph i kiểm định, ph ơng pháp tính toán hệ số t ơng quan Kendall đòi hỏi th i gian của máy tính khá cao. Tuy nhiên, nếu một dữ liệu d ới 5000 đối t ợng thì một máy vi tính có thể tính toán khá dễ dàng. R dùng hàm cor.test với thông số method=”kendall” để ớc tính hệ số t ơng quan Kendall: > cor.test(age, chol, method="kendall") Kendall's rank correlation tau data: age and chol z = 4.755, p-value = 1.984e-06 alternative hypothesis: true tau is not equal to 0 sample estimates: tau 0.8333333 74 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Warning message: Cannot compute exact p-value with ties in: cor.test.default(age, chol, method = "kendall") 10.2 Mô hình của hồi qui tuyến tính đơn giản Để tiện việc theo dõi và mô t mô hình, gọi độ tuổi cho cá nhân i là xi và cholesterol là yi. đây i = 1, 2, 3, …, 18. Mô hình hồi tuyến tính phát biểu rằng: yi = α + β xi + ε i Nói cách khác, ph ơng trình trên gi định rằng độ cholesterol của một cá nhân bằng một hằng số α cộng với một hệ số β liên quan đến độ tuổi, và một sai số εi. Trong ph ơng trình trên, α là chặn (intercept, tức giá trị lúc xi =0), và β là độ dốc (slope hay gradient). Trong thực tế, α và β là hai thông số (paramater, còn gọi là regression coefficient hay hệ số hồi qui), và εi là một biến số theo luật phân phối chuẩn với trung bình 0 và ph ơng sai σ2 . Các thông số α, β và σ2 ph i đ ợc ớc tính từ dữ liệu. Ph ơng pháp để ớc tính các thông số này là ph ơng pháp bình phương nhỏ nhất (least squares method). Nh tên gọi, ph ơng pháp bình ph ơng nhỏ nhất tìm giá trị α, β sao cho ∑  y − (α + β x ) n i =1 i i 2 nhỏ nhất. Sau vài thao tác toán, có thể chứng minh dễ dàng rằng, ớc số cho α và β đáp ứng điều kiện đó là: ∑ ( x − x )( y − y ) n βˆ = i =1 ∑(x − x ) i i n i =1 2 ) ) và α = y − β x i ) ) đây, x và y là giá trị trung bình của biến số x và y. Chú ý, tôi viết α và β (với dấu mũ phía trên) là để nhắc nh rằng đây là hai ớc số (estimates) của α và β, chứ không ph i α và β (chúng ta không biết chính xác α và β, nh ng chỉ có thể ớc tính mà thôi). ) ) Sau khi đã có ớc số α và β , chúng ta có thể ớc tính độ cholesterol trung bình cho từng độ tuổi nh sau: ) yi = αˆ + βˆ xi Tất nhiên, yˆi đây chỉ là số trung bình cho độ tuổi xi, và phần còn l i (tức yi - yˆi ) gọi là phần dư (residual). Và ph ơng sai của phần d có thể ớc tính nh sau: ∑ ( y − yˆ ) n s2 = i =1 i n−2 i . đây, s2 chính là ớc số của σ2. 75 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn ) Hàm lm (viết tắt từ linear model) trong R có thể tính toán các giá trị của α và β , cũng nh s2 một cách nhanh gọn. Chúng ta tiếp tục với ví dụ bằng R nh sau: ) > lm(chol ~ age) Call: lm(formula = chol ~ age) Coefficients: (Intercept) 1.08922 age 0.05779 Trong lệnh trên, “chol ~ age” có nghĩa là mô t chol là một hàm số của age. Kết ) ) qu tính toán của lm cho thấy α = 1.0892 và β = 0.05779. Nói cách khác, với hai thông số này, chúng ta có thể ớc tính độ cholesterol cho bất cứ độ tuổi nào trong kho ng tuổi của mẫu bằng ph ơng trình tuyến tính: yˆi = 1.08922 + 0.05779 x age Ph ơng trình này có nghĩa là khi độ tuổi tăng 1 năm thì độ cholesterol tăng kho ng 0.058 mmol/L. Thật ra, hàm lm còn cung cấp cho chúng ta nhiều thông tin khác, nh ng chúng ta ph i đ a các thông tin này vào một object. Gọi object đó là reg, thì lệnh sẽ là: > reg <- lm(chol ~ age) > summary(reg) Call: lm(formula = chol ~ age) Residuals: Min 1Q Median -0.40729 -0.24133 -0.04522 3Q 0.17939 Max 0.63040 Coefficients: Estimate Std. Error t value (Intercept) 1.089218 0.221466 4.918 age 0.057788 0.005399 10.704 --Signif. codes: 0 '***' 0.001 '**' 0.01 Pr(>|t|) 0.000154 *** 1.06e-08 *** '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3027 on 16 degrees of freedom Multiple R-Squared: 0.8775, Adjusted R-squared: 0.8698 F-statistic: 114.6 on 1 and 16 DF, p-value: 1.058e-08 Lệnh thứ hai, summary(reg), yêu cầu R liệt kê các thông tin tính toán trong reg. Phần kết qu chia làm 3 phần: 76 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn (a) Phần 1 mô t phần d (residuals) của mô hình hồi qui: Residuals: Min 1Q Median -0.40729 -0.24133 -0.04522 3Q 0.17939 Max 0.63040 Chúng ta biết rằng trung bình phần d ph i là 0, và đây, số trung vị là -0.04, cũng không xa 0 bao nhiêu. Các số quantiles 25% (1Q) và 75% (3Q) cũng khá cân đối chung quan số trung vị, cho thấy phần d của ph ơng trình này t ơng đối cân đối. ) ) (b) Phần hai trình bày ớc số của α và β cùng với sai số chuẩn và giá trị của kiểm định t. ) Giá trị kiểm định t cho β là 10.74 với trị số p = 1.06e-08, cho thấy β không ph i bằng 0. Nói cách khác, chúng ta có bằng chứng để cho rằng có một mối liên hệ giữa cholesterol và độ tuổi, và mối liên hệ này có ý nghĩa thống kê. Coefficients: Estimate Std. Error t value (Intercept) 1.089218 0.221466 4.918 age 0.057788 0.005399 10.704 --Signif. codes: 0 '***' 0.001 '**' 0.01 Pr(>|t|) 0.000154 *** 1.06e-08 *** '*' 0.05 '.' 0.1 ' ' 1 (c) Phần ba của kết qu cho chúng ta thông tin về ph ơng sai của phần d (residual mean square). đây, s2 = 0.3027. Trong kết qu này còn có kiểm định F, cũng chỉ là một kiểm định xem có qu thật β bằng 0, tức có ý nghĩa t ơng tự nh kiểm định t trong phần trên. Nói chung, trong tr ng hợp phân tích hồi qui tuyến tính đơn gi n (với một yếu tố) chúng ta không cần ph i quan tâm đến kiểm định F. Residual standard error: 0.3027 on 16 degrees of freedom Multiple R-Squared: 0.8775, Adjusted R-squared: 0.8698 F-statistic: 114.6 on 1 and 16 DF, p-value: 1.058e-08 Ngoài ra, phần 3 còn cho chúng ta một thông tin quan trọng, đó là trị số R2 hay hệ số xác định bội (coefficient of determination). Tức là bằng tổng bình ph ơng giữa số ớc tính và trung bình chia cho tổng bình ph ơng số quan sát và trung bình. Trị số R2 trong ví dụ này là 0.8775, có nghĩa là ph ơng trình tuyến tính (với độ tuổi là một yếu tố) gi i thích kho ng 88% các khác biệt về độ cholesterol giữa các cá nhân. Tất nhiên trị số R2 có giá trị từ 0 đến 100% (hay 1). Giá trị R2 càng cao là một dấu hiệu cho thấy mối liên hệ giữa hai biến số độ tuổi và cholesterol càng chặt chẽ. Một hệ số cũng cần đề cập đây là hệ số điều chỉnh xác định bội (mà trong kết qu trên R gọi là “Adjusted R-squared”). Đây là hệ số cho chúng ta biết mức độ c i tiến của ph ơng sai phần d (residual variance) do yếu tố độ tuổi có mặt trong mô hình tuyến tính. Nói chung, hệ số này không khác mấy so với hệ số xác định bội, và chúng ta cũng không cần chú tâm quá mức. Giả định của phân tích hồi qui tuyến tính 77 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Tất c các phân tích trên dựa vào một số gi định quan trọng nh sau: (a) x là một biến số cố định hay fixed, (“cố định” nhiên trong đo l ng); đây có nghĩa là không có sai sót ngẫu (b) εi phân phối theo luật phân phối chuẩn; (c) εi có giá trị trung bình (mean) là 0; (d) εi có ph ơng sai σ2 cố định cho tất c xi; và (e) các giá trị liên tục của εi không có liên hệ t ơng quan với nhau (nói cách khác, ε1 và ε2 không có liên hệ với nhau). Nếu các gi định này không đ ợc đáp ứng thì ph ơng trình mà chúng ta ớc tính có vấn đề hợp lí (validity). Do đó, tr ớc khi trình bày và diễn dịch mô hình trên, chúng ta cần ph i kiểm tra xem các gi định trên có đáp ứng đ ợc hay không. Trong tr ng hợp này, gi định (a) không ph i là vấn đề, vì độ tuổi không ph i là một biến số ngẫu nhiên, và không có sai số khi tính độ tuổi của một cá nhân. Đối với các gi định (b) đến (e), cách kiểm tra đơn gi n nh ng hữu hiệu nhất là bằng cách xem xét mối liên hệ giữa yˆi , xi , và phần d ei ( ei = yi − yˆi ) bằng những đồ thị tán x . Với lệnh fitted() chúng ta có thể tính toán yˆi cho từng cá nhân nh sau (ví dụ đối với cá nhân 1, 46 tuổi, độ cholestrol có thể tiên đoán nh sau: 1.08922 + 0.05779 x 46 = 3.747). > fitted(reg) 1 2 3 4 5 6 7 8 3.747483 2.244985 4.094214 2.822869 4.383156 2.533927 2.707292 3.169600 9 10 11 12 13 14 15 16 2.360562 3.574118 4.383156 2.996234 2.360562 4.729886 3.400753 3.863060 17 18 2.707292 3.920849 Với lệnh resid() chúng ta có thể tính toán phần d ei cho từng cá nhân nh sau (với đối t ợng 1, e1 = 3.5 – 3.74748 = -0.24748): > resid(reg) 1 2 3 4 5 -0.247483426 -0.344985415 -0.094213736 -0.222869265 0.116844338 7 8 9 10 11 0.192707505 0.630400424 -0.260562185 0.225881729 -0.283155662 13 14 15 16 17 0.139437815 -0.129885972 -0.200753116 0.336939804 -0.407292495 78 6 0.466072660 12 0.003765579 18 0.079151419 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Để kiểm tra các gi định trên, chúng ta có thể vẽ một lo t 4 đồ thị mà tôi sẽ gi i thích sau đây: #yêu cầu R dành ra 4 cửa sổ #vẽ các đồ thị trong reg > op <- par(mfrow=c(2,2)) > plot(reg) Normal Q-Q 8 1 0.2 -1 0.0 6 0 Standardized residuals 2 8 6 -0.4 Residuals 0.4 0.6 Residuals vs Fitted 17 17 3.0 3.5 4.0 4.5 -2 -1 0 1 2 Fitted values Theoretical Quantiles Scale-Location Residuals vs Leverage 1 8 2 8 0.5 1 6 0 1.0 17 -1 Standardized residuals 6 0.5 Standardized residuals 1.5 2.5 2 0.0 Cook's distance 2.5 3.0 3.5 4.0 4.5 0.00 0.05 Fitted values 0.10 0.5 0.15 0.20 0.25 Leverage Biểu đồ 19. Phân tích phần d để kiểm tra các gi định trong phân tích hồi qui tuyến tính. (a) Đồ thị bên trái dòng 1 vẽ phần d ei và giá trị tiên đoán cholesterol yˆi . Đồ thị này cho thấy các giá trị phần d tập chung quanh đ ng y = 0, cho nên gi định (c), hay εi có giá trị trung bình 0, là có thể chấp nhận đ ợc. (b) Đồ thị bên ph i dòng 1 vẽ giá trị phần d và giá trị kì vọng dựa vào phân phối chuẩn. Chúng ta thấy các số phần d tập trung rất gần các giá trị trên đ ng chuẩn, và do đó, gi định (b), tức εi phân phối theo luật phân phối chuẩn, cũng có thể đáp ứng. (c) Đồ thị bên trái dòng 2 vẽ căn số phần d chuẩn (standardized residual) và giá trị của yˆi . Đồ thị này cho thấy không có gì khác nhau giữa các số phần d chuẩn cho các giá trị 79 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn của yˆi , và do đó, gi định (d), tức εi có ph ơng sai σ2 cố định cho tất c xi, cũng có thể đáp ứng. Nói chung qua phân tích phần d , chúng ta có thể kết luận rằng mô hình hồi qui tuyến tính mô t mối liên hệ giữa độ tuổi và cholesterol một cách khá đầy đủ và hợp lí. Mô hình tiên đoán Sau khi mô hình tiên đoán cholesterol đã đ ợc kiểm tra và tính hợp lí đã đ ợc thiết lập, chúng ta có thể vẽ đ ng biểu diễn của mối liên hệ giữa độ tuổi và cholesterol bằng lệnh abline nh sau (xin nhắc l i object của phân tích là reg): 2.0 2.5 3.0 chol 3.5 4.0 4.5 > plot(chol ~ age, pch=16) > abline(reg) 20 30 40 50 60 age Biểu đồ 20. Đ cholesterol. ng biểu diễn mối liên hệ giữa độ tuổi (age) và ) ) Nh ng mỗi giá trị yˆi đ ợc tính từ ớc số α và β , mà các ớc số này đều có sai số chuẩn, cho nên giá trị tiên đoán yˆi cũng có sai số. Nói cách khác, yˆi chỉ là trung bình, 80 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn nh ng trong thực tế có thể cao hơn hay thấp hơn tùy theo chọn mẫu. Kho ng tin cậy 95% này có thể ớc tính qua R bằng các lệnh sau đây: reg <- lm(chol ~ age) new <- data.frame(age = seq(15, 70, 5)) pred.w.plim <- predict.lm(reg, new, interval="prediction") pred.w.clim <- predict.lm(reg, new, interval="confidence") resc <- cbind(pred.w.clim, new) resp <- cbind(pred.w.plim, new) plot(chol ~ age, pch=16) lines(resc$fit ~ resc$age) lines(resc$lwr ~ resc$age, col=2) lines(resc$upr ~ resc$age, col=2) lines(resp$lwr ~ resp$age, col=4) lines(resp$upr ~ resp$age, col=4) 2.0 2.5 3.0 chol 3.5 4.0 4.5 > > > > > > > > > > > > 20 30 40 50 60 age Biểu đồ 21. Giá trị tiên đoán và kho ng tin cậy 95%. Biểu đồ trên vẽ giá trị tiên đoán trung bình yˆi (đ ng thẳng màu đen), và kho ng tin cậy 95% của giá trị này là đ ng màu đỏ. Ngoài ra, đ ng màu xanh là kho ng tin cậy của giá trị tiên đoán cholesterol cho một độ tuổi mới trong quần thể. 81 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 10.3 Mô hình hồi qui tuyến tính đa biến (multiple linear regression) Mô hình đ ợc diễn đ t qua ph ơng trình yi = α + β xi + ε i có một yếu tố duy nhất (đó là x), và vì thế th ng đ ợc gọi là mô hình hồi qui tuyến tính đơn gi n (simple linear regression model). Trong thực tế, chúng ta có thể phát triển mô hình này thành nhiều biến, chứ không chỉ giới h n một biến nh trên, chẳng h n nh : yi = α + β1 x1i + β 2 x2i + ... + β k xki + ε i Chú ý trong ph ơng trình trên, chúng ta có nhiều biến x (x1, x2, … đến xk), và mỗi biến có một thông số β j (j = 1, 2, …, k) cần ph i ớc tính. Vì thế mô hình này còn đ ợc gọi là mô hình hồi qui tuyến tính đa biến. Ví dụ 16. Chúng ta quay l i nghiên cứu về mối liên hệ giữa độ tuổi, bmi và cholesterol. Trong ví dụ, chúng ta chỉ mới xét mối liên hệ giữa độ tuổi và cholesterol, mà ch a xem đến mối liên hệ giữa c hai yếu tố độ tuổi và bmi và cholesterol. Biểu đồ sau đây cho chúng ta thấy mối liên hệ giữa ba biến số này: > pairs(data) 22 24 26 50 60 20 24 26 20 30 40 age chol 20 30 40 50 60 2.0 2.5 3.0 3.5 4.0 4.5 20 22 bmi 2.0 2.5 3.0 3.5 4.0 4.5 Biểu đồ 22. Giá trị tiên đoán và kho ng tin cậy 95%. Cũng nh giữa độ tuổi và cholesterol, mối liên hệ giữa bmi và cholesterol cũng gần tuân theo một đ ng thằng. Biểu đồ trên còn cho chúng ta thấy độ tuổi và bmi có liên hệ với 82 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn nhau. Thật vậy, phân tích hồi qui tuyến tính đơn gi n giữa bmi và cholesterol cho thấy nh mối liên hệ này có ý nghĩa thống kê: > summary(lm(chol ~ bmi)) Call: lm(formula = chol ~ bmi) Residuals: Min 1Q Median -0.9403 -0.3565 -0.1376 3Q 0.3040 Max 1.4330 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.83187 1.60841 -1.761 0.09739 . bmi 0.26410 0.06861 3.849 0.00142 ** --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.623 on 16 degrees of freedom Multiple R-Squared: 0.4808, Adjusted R-squared: 0.4483 F-statistic: 14.82 on 1 and 16 DF, p-value: 0.001418 BMI gi i thích kho ng 48% độ dao động về cholesterol giữa các cá nhân. Nh ng vì BMI cũng có liên hệ với độ tuổi, chúng ta muốn biết nếu hai yếu tố này đ ợc phân tích cùng một lúc thì yếu tố nào quan trọng hơn. Để biết nh h ng của c hai yếu tố age (x1) và bmi (t m gọi là x2) đến cholesterol (y) qua một mô hình hồi qui tuyến tính đa biến, và mô hình đó là: yi = α + β1 x1i + β 2 x2i + ε i hay ph ơng trình cũng có thể mô t bằng kí hiệu ma trận: Y = Xβ + ε mà tôi vừa trình bày trên. đây, Y là một vector vector 18 x 1, X là một matrix 18 x 2 phần tử, β và một vector 2 x 1, và ε là vector gồm 18 x 1 phần tử. Để ớc tính hai hệ số hồi qui, β1 và β2 chúng ta cũng ứng dụng hàm lm() trong R nh sau: > mreg <- lm(chol ~ age + bmi) > summary(mreg) Call: lm(formula = chol ~ age + bmi) Residuals: Min 1Q Median -0.3762 -0.2259 -0.0534 3Q 0.1698 Max 0.5679 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.455458 0.918230 0.496 0.627 83 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn age 0.054052 0.007591 7.120 3.50e-06 *** bmi 0.033364 0.046866 0.712 0.487 --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3074 on 15 degrees of freedom Multiple R-Squared: 0.8815, Adjusted R-squared: 0.8657 F-statistic: 55.77 on 2 and 15 DF, p-value: 1.132e-07 Kết qu phân tích trên cho thấy ớc số α̂ = 0.455, β̂1 = 0.054 và β̂ 2 = 0.0333. Nói cách khác, chúng ta có ph ơng trình ớc đoán độ cholesterol dựa vào hai biến số độ tuổi và bmi nh sau: Cholesterol = 0.455 + 0.054(age) + 0.0333(bmi) Ph ơng trình cho biết khi độ tuổi tăng 1 năm thì cholesterol tăng 0.054 mg/L ( ớc số này không khác mấy so với 0.0578 trong ph ơng trình chỉ có độ tuổi), và mỗi 1 kg/m2 tăng BMI thì cholesterol tăng 0.0333 mg/L. Hai yếu tố này “gi i thích” kho ng 88.2% (R2 = 0.8815) độ dao động của cholesterol giữa các cá nhân. Chúng ta chú ý ph ơng trình với độ tuổi (trong phân tích phần tr ớc) gi i thích kho ng 87.7% độ dao động cholesterol giữa các cá nhân. Khi chúng ta thêm yếu tố BMI, hệ số này tăng lên 88.2%, tức chỉ 0.5%. Câu hỏi đặt ra là 0.5% tăng tr ng này có ý nghĩa thống kê hay không. Câu tr l i có thể xem qua kết qu kiểm định yếu tố bmi với trị số p = 0.487. Nh vậy, bmi không cung cấp cho chúng thêm thông tin hay tiên đoán cholesterol hơn những gì chúng ta đã có từ độ tuổi. Nói cách khác, khi độ tuổi đã đ ợc xem xét, thì nh h ng của bmi không còn ý nghĩa thống kê. Điều này có thể hiểu đ ợc, b i vì qua Biểu đồ 10.5 chúng ta thấy độ tuổi và bmi có một mối liên hệ khá cao. Vì hai biến này có t ơng quan với nhau, chúng ta không cần c hai trong ph ơng trình. (Tuy nhiên, ví dụ này chỉ có tính cách minh họa cho việc tiến hành phân tích hồi qui tuyến tính đa biến bằng R, chứ không có ý định mô phỏng dữ liệu theo định h ớng sinh học). 84 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 3.0 4.0 2.0 1.0 -2 -1 0 1 2 Residuals vs Leverage 3.5 4.0 0.5 1 16 0 0.8 0.4 3.0 8 -1 16 2 Scale-Location Standardized residuals Theoretical Quantiles 6 2.5 0.0 4.5 Fitted values 8 1.2 3.5 0.0 Standardized residuals 2.5 8 16 6 -1.0 0.0 0.4 16 -0.4 Residuals 8 6 Normal Q-Q Standardized residuals Residuals vs Fitted 4.5 Cook's distance15 0.00 0.10 Fitted values 0.20 0.30 Leverage Biểu đồ 23. Phân tích phần d để kiểm tra các gi định trong phân tích hồi qui tuyến tính đa biến. Tuy BMI không có ý nghĩa thống kê trong tr ng hợp này, Biểu đồ 10.6 cho thấy các gi định về mô hình hồi qui tuyến tính có thể đáp ứng. 11. Phân tích ph ơng sai 11.1 Phân tích ph ơng sai đơn giản (one-way analysis of variance ANOVA) Ví dụ 17. B ng d ới đây so sánh độ galactose trong 3 nhóm bệnh nhân: nhóm 1 gồm 9 bệnh nhân với bệnh Crohn; nhóm 2 gồm 11 bệnh nhân với bệnh viêm ruột kết (colitis); và nhóm 3 gồm 20 đối t ợng không có bệnh (gọi là nhóm đối chứng). Câu hỏi đặt ra là độ galactose giữa 3 nhóm bệnh nhân có khác nhau hay không? Độ galactose cho 3 nhóm bệnh nhân Crohn, viêm ruột kết và đối ch ng Nhóm 1: bệnh Crohn Nhóm 2: bệnh viêm ruột kết Nhóm 3: đối chứng (control) 85 Phân tích số liệu và biểu đồ bằng R 1343 1393 1420 1641 1897 2160 2169 2279 2890 Nguyễn Văn Tuấn 1264 1314 1399 1605 2385 2511 2514 2767 2827 2895 3011 1809 1926 2283 2384 2447 2479 2495 2525 2541 2769 2850 2964 2973 3171 3257 3271 3288 3358 3643 3657 n=9 n=11 n=20 Trung bình: 1910 Trung bình: 2226 Trung bình: 2804 SD: 516 SD: 727 SD: 527 Chú thích: SD là độ lệch chuẩn (standard deviation). Gọi giá trị trung bình của ba nhóm là µ1, µ2, và µ3, và nói theo ngôn ngữ của kiểm định gi thiết thì gi thiết đ o là: Và gi thiết chính là: Ho: µ1 = µ2 = µ3 HA: có một khác biệt giữa 3 µj (j = 1,2,3) Tho t đầu có lẽ b n đọc, sau khi đã học qua ph ơng pháp so sánh hai nhóm bằng kiểm định t, sẽ nghĩ rằng chúng ta cần làm 3 so sánh bằng kiểm định t: giữa nhóm 1 và 2, nhóm 2 và 3, và nhóm 1 và 3. Nh ng ph ơng pháp này không hợp lí, vì có ba ph ơng sai khác nhau. Ph ơng pháp thích hợp cho so sánh là phân tích ph ơng sai. Phân tích ph ơng sai có thể ứng dụng để so sánh nhiều nhóm cùng một lúc (simultaneous comparisons). Để minh họa cho ph ơng pháp phân tích ph ơng sai, chúng ta ph i dùng kí hiệu. Gọi độ galactose của bệnh nhân i thuộc nhóm j (j = 1, 2, 3) là xij. Mô hình phân tích ph ơng sai phát biểu rằng: xij = µ + α i + ε ij Hay cụ thể hơn: xi1 = µ + α1 + εi1 xi2 = µ + α2 + εi2 xi3 = µ + α3 + εi3 Tr ớc hết, chúng ta cần ph i nhập dữ liệu vào R. B ớc thứ nhất là báo cho R biết rằng chúng ta có ba nhóm bệnh nhân (1, 2 v ), nhóm 1 gồm 9 ng i, nhóm 2 có 11 ng i, và nhóm 3 có 20 ng i: > group <- c(1,1,1,1,1,1,1,1,1, 2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3) 86 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Để phân tích ph ơng sai, chúng ta ph i định nghĩa biến group là một yếu tố - factor. > group <- as.factor(group) B ớc kế tiếp, chúng ta n p số liệu galactose cho từng nhóm nh định nghĩa trên (gọi object là galactose): > galactose <- c(1343,1393,1420,1641,1897,2160,2169,2279,2890, 1264,1314,1399,1605,2385,2511,2514,2767,2827,2895,3011, 1809,2850,1926,2964,2283,2973,2384,3171,2447,3257,2479,3271,2495,3288, 2525,3358,2541,3643,2769,3657) Đ a hai biến group và galactose vào một dataframe và gọi là data: > data <- data.frame(group, galactose) > attach(data) Sau khi đã có dử liệu sẵn sàng, chúng ta dùng hàm lm() để phân tích ph ơng sai nh sau: > analysis <- lm(galactose ~ group) Trong hàm trên chúng ta cho R biết biến galactose là một hàm số của group. Gọi kết qu phân tích là analysis. Kết quả phân tích phương sai. Bây gi chúng ta dùng lệnh anova để biết kết qu phân tích: > anova(analysis) Analysis of Variance Table Response: galactose Df Sum Sq Mean Sq F value Pr(>F) group 2 5683620 2841810 8.6655 0.0008191 *** Residuals 37 12133923 327944 --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Trong kết qu trên, có ba cột: Df (degrees of freedom) là bậc tự do; Sum Sq là tổng bình ph ơng (sum of squares), Mean Sq là trung bình bình ph ơng (mean square); F value là giá trị F; và Pr(>F) là trị số P liên quan đến kiểm định F. 11.2 So sánh nhiều nhóm (multiple comparisons) và điều chỉnh trị số p Cho k nhóm, chúng ta có ít nhất là k(k-1)/2 so sánh. Ví dụ trên có 3 nhóm, cho nên tổng số so sánh kh dĩ là 3 (giữa nhóm 1 và 2, nhóm 1 và 3, và nhóm 2 và 3). Khi k=10, số lần so sánh có thể lên rất cao. Nh đã đề cập trong ch ơng 7, khi có nhiều so sánh, trị số p tính toán từ các kiểm định thống kê không còn ý nghĩa ban đầu nữa, b i vì các kiểm định này có thể cho ra kết qu d ơng tính gi (tức kết qu với p<0.05 nh ng 87 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn trong thực tế không có khác nhau hay nh h ng). Do đó, trong tr sánh, chúng ta cần ph i điều chỉnh trị số p sao cho hợp lí. ng hợp có nhiều so Có khá nhiều ph ơng pháp điều chỉnh trị số p, và 4 ph ơng pháp thông dụng nhất là: Bonferroni, Scheffé, Holm và Tukey (tên của 4 nhà thống kê học danh tiếng). Ph ơng pháp nào thích hợp nhất? Không có câu tr l i dứt khoát cho câu hỏi này, nh ng hai điểm sau đây có thể giúp b n đọc quyết định tốt hơn: (a) Nếu k < 10, chúng ta có thể áp dụng bất cứ ph ơng pháp nào để điều chỉnh trị số p. Riêng cá nhân tôi thì thấy ph ơng pháp Tukey th ng rất hữu ích trong so sánh. (b) Nếu k>10, ph ơng pháp Bonferroni có thể tr nên rất “b o thủ”. B o thủ đây có nghĩa là ph ơng pháp này rất ít khi nào tuyên bố một so sánh có ý nghĩa thống kê, dù trong thực tế là có thật! Trong tr ng hợp này, hai ph ơng pháp Tukey, Holm và Scheffé có thể áp dụng. Quay l i ví dụ trên, các trị số p trên đây là những trị số ch a đ ợc điều chỉnh cho so sánh nhiều lần. Trong ch ơng về trị số p, tôi đã nói các trị số này phóng đ i ý nghĩa thống kê, không ph n ánh trị số p lúc ban đầu (tức 0.05). Để điều chỉnh cho nhiều so sánh, chúng ta ph i sử dụng đến ph ơng pháp điều chỉnh Bonferroni. Chúng ta có thể dùng lệnh pairwise.t.test để có đ ợc tất c các trị số p so sánh giữa ba nhóm nh sau: > pairwise.t.test(galactose, group, p.adj="bonferroni") Pairwise comparisons using t tests with pooled SD data: galactose and group 1 2 2 0.6805 3 0.0012 0.0321 P value adjustment method: bonferroni Kết qu trên cho thấy trị số p giữa nhóm 1 (Crohn) và viêm ruột kết là 0.6805 (tức không có ý nghĩa thống kê); giữa nhóm Crohn và đối chứng là 0.0012 (có ý nghĩa thống kê), và giữa nhóm viêm ruột kết và đối chứng là 0.0321 (tức cũng có ý nghĩa thống kê). Một ph ơng pháp điều chỉnh trị số p khác có tên là ph ơng pháp Holm: > pairwise.t.test(galactose, group) Pairwise comparisons using t tests with pooled SD data: galactose and group 88 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 1 2 2 0.2268 3 0.0012 0.0214 P value adjustment method: holm Kết qu này cũng không khác so với ph ơng pháp Bonferroni. Tất c các ph ơng pháp so sánh trên sử dụng một sai số chuẩn chung cho c ba nhóm. Nếu chúng ta muốn sử dụng cho từng nhóm thì lệnh sau đây (pool.sd=F) sẽ đáp ứng yêu cầu đó: > pairwise.t.test(galactose, group, pool.sd=FALSE) Pairwise comparisons using t tests with non-pooled SD data: galactose and group 1 2 2 0.2557 3 0.0017 0.0544 P value adjustment method: holm Một lần nữa, kết qu này cũng không làm thay đổi kết luận. Trong các ph ơng pháp trên, chúng ta chỉ biết trị số p so sánh giữa các nhóm, nh ng không biết mức độ khác biệt cũng nh kho ng tin cậy 95% giữa các nhóm. Để có những ớc số này, chúng ta cần đến một hàm khác có tên là aov (viết tắt từ analysis of variance) và hàm TukeyHSD (HSD là viết tắt từ Honest Significant Difference, t m dịch nôm na là “Khác biệt có ý nghĩa thành thật”) nh sau: > res <- aov(galactose ~ group) > TukeyHSD (res) Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = galactose ~ group) $group diff lwr upr p adj 2-1 316.3232 -312.09857 944.745 0.4439821 3-1 894.2778 333.07916 1455.476 0.0011445 3-2 577.9545 53.11886 1102.790 0.0281768 Kết qu trên cho chúng ta thấy nhóm 3 và 1 khác nhau kho ng 894 đơn vị, và kho ng tin cậy 95% từ 333 đến 1455 đơn vị. T ơng tự, galactose trong nhóm bệnh nhân viêm ruột kết thấp hơn nhóm đối chứng (nhóm 3) kho ng 578 đơn vị, và kho ng tin cậy 95% từ 53 đến 1103. 89 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 3-2 3-1 2-1 95% family-wise confidence level 0 500 1000 1500 Differences in mean levels of group Biểu đồ 24. Trung bình hiệu và kho ng tin cậy 95% giữa nhóm 1 và 2, 1 và 3, và 3 và 2. Trục hoành là độ galactose, trục tung là ba so sánh. 11.3 Phân tích bằng ph ơng pháp phi tham số Ph ơng pháp so sánh nhiều nhóm phi tham số (non-parametric statistics) t ơng đ ơng với ph ơng pháp phân tích ph ơng sai là Kruskal-Wallis. Cũng nh ph ơng pháp Wilcoxon so sánh hai nhóm theo ph ơng pháp phi tham số, ph ơng pháp Kruskal-Wallis cũng biến đổi số liệu thành thứ bậc (ranks) và phân tích độ khác biệt thứ bậc này giữa các nhóm. Hàm kruskal.test trong R có thể giúp chúng ta trong kiểm định này: > kruskal.test(galactose ~ group) Kruskal-Wallis rank sum test data: galactose by group Kruskal-Wallis chi-squared = 12.1381, df = 2, p-value = 0.002313 Trị số p từ kiểm định này khá thấp (p = 0.002313) cho thấy có sự khác biệt giữa ba nhóm nh phân tích ph ơng sai qua hàm lm trên đây. Tuy nhiên, một bất tiện của kiểm định phi tham số Kruskal-Wallis là ph ơng pháp này không cho chúng ta biết hai nhóm nào khác nhau, mà chỉ cho một trị số p chung. Trong nhiều tr ng hợp, phân tích phi tham số nh kiểm định Kruskal-Wallis th ng không có hiệu qu nh các ph ơng pháp thống kê tham số (parametric statistics). 90 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 11.4 Phân tích ph ơng sai hai chiều (two-way analysis of variance ANOVA) Phân tích ph ơng sai đơn gi n hay một chiều chỉ có một yếu tố (factor). Nh ng phân tích ph ơng sai hai chiều (two-way ANOVA), nh tên gọi, có hai yếu tố. Ph ơng pháp phân tích ph ơng sai hai chiều chỉ đơn gi n khai triển từ ph ơng pháp phân tích ph ơng sai đơn gi n. Thay vì ớc tính ph ơng sai của một yếu tố, ph ơng pháp phân sai hai chiều ớc tính ph ơng sai của hai yếu tố. Ví dụ 18. Trong ví dụ sau đây, để đánh giá hiệu qu của một kĩ thuật sơn mới, các nhà nghiên cứu áp dụng sơn trên 3 lo i vật liệu (1, 2 v 3) trong hai điều kiện (1, 2). Mỗi điều kiện và lo i vật liệu, nghiên cứu đ ợc lặp l i 3 lần. Độ bền đ ợc đo là chỉ số bền bĩ (t m gọi là score). Tổng cộng, có 18 số liệu nh sau: Độ bền bĩ c a sơn cho 2 điều kiện và 3 vật liệu Điều kiện (i) 1 2 1 4.1, 3.9, 4.3 2.7, 3.1, 2.6 Vật liệu (j) 2 3.1, 2.8, 3.3 1.9, 2.2, 2.3 3 3.5, 3.2, 3.6 2.7, 2.3, 2.5 Gọi xij là score của điều kiện i (i = 1, 2) cho vật liệu j (j = 1, 2, 3). (Để đơn gi n hóa vấn đề, chúng ta t m th i bỏ qua k đối t ợng). Mô hình phân tích ph ơng sai hai chiều phát biểu rằng: xij = µ + α i + β j + ε ij µ là số trung bình cho toàn quần thể, các hệ số αi ( nh h ng của điều kiện i)và βj ( nh h ng của vật liệu j) cần ph i ớc tính từ số liệu thực tế. εij đ ợc gi định tuân theo luật phân phối chuẩn với trung bình 0 và ph ơng sai σ2. Để phân tích bằng R, chúng ta cần ph i tổ chức dữ liệu sao cho có 4 biến nh sau: Condition (điều kiện) 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 Material (vật liệu) 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 Đối tượng Score 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 4.1 3.9 4.3 3.1 2.8 3.3 3.5 3.2 3.6 2.7 3.1 2.6 1.9 2.2 2.3 91 Phân tích số liệu và biểu đồ bằng R 2 2 2 3 3 3 Nguyễn Văn Tuấn 16 17 18 2.7 2.3 2.5 Chúng ta có thể t o ra một dãy số bằng cách sử dụng hàm gl (generating levels). > condition <- gl(2, 9, 18) > material <- gl(3, 3, 18) Và t o nên 18 mã số (tứ 1 đến 18): > id <- 1:18 Sau cùng là số liệu cho score: > score <- c(4.1,3.9,4.3, 3.1,2.8,3.3, 3.5,3.2,3.6, 2.7,3.1,2.6, 1.9,2.2,2.3, 2.7,2.3,2.5) Tất c cho vào một dataframe tên là data: > data <- data.frame(condition, material, id, score) > attach(data) Bây gi số liệu đã sẵn sàng cho phân tích. Để phân tích ph ơng sai hai chiều, chúng ta vẫn sử dụng lệnh lm với các thông số nh sau: > twoway <- lm(score ~ condition + material) > anova(twoway) Analysis of Variance Table Response: score Df Sum Sq Mean Sq F value Pr(>F) condition 1 5.0139 5.0139 95.575 1.235e-07 *** material 2 2.1811 1.0906 20.788 6.437e-05 *** Residuals 14 0.7344 0.0525 --Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Ba nguồn dao động (variation) của score đ ợc phân tích trong b ng trên. Qua trung bình bình ph ơng (mean square), chúng ta thấy nh h ng của điều kiện có vẻ quan trọng hơn là nh h ng của vật liệu thí nghiệm. Tuy nhiên, c hai nh h ng đều có ý nghĩa thống kê, vì trị số p rất thấp cho hai yếu tố. Chúng ta yêu cầu R tóm l ợc các ớc số phân tích bằng lệnh summary: > summary(twoway) Call: lm(formula = score ~ condition + material) Residuals: Min 1Q -0.32778 -0.16389 Median 0.03333 3Q 0.16111 Max 0.32222 92 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Coefficients: Estimate Std. Error t value (Intercept) 3.9778 0.1080 36.841 condition2 -1.0556 0.1080 -9.776 material2 -0.8500 0.1322 -6.428 material3 -0.4833 0.1322 -3.655 --Signif. codes: 0 '***' 0.001 '**' 0.01 Pr(>|t|) 2.43e-15 1.24e-07 1.58e-05 0.0026 *** *** *** ** '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.229 on 14 degrees of freedom Multiple R-Squared: 0.9074, Adjusted R-squared: 0.8875 F-statistic: 45.72 on 3 and 14 DF, p-value: 1.761e-07 Kết qu trên cho thấy so với điều kiện 1, điều kiện 2 có score thấp hơn kho ng 1.056 và sai số chuẩn là 0.108, với trị số p = 1.24e-07, tức có ý nghĩa thống kê. Ngoài ra, so với vật liệu 1, score cho vật liệu 2 và 3 cũng thấp hơn đáng kể với độ thấp nhất ghi nhận vật liệu 2, và nh h ng của vật liệu thí nghiệm cũng có ý nghĩa thống kê. Giá trị có tên là “Residual standard error” đ ợc ớc tính từ trung bình bình ph ơng phần d trong phần (a), tức là 0.0525 = 0.229, tức là ớc số của σˆ . Hệ số xác định bội (R2) cho biết hai yếu tố điều kiện và vật liệu gi i thích kho ng 91% độ dao động của toàn bộ mẫu. Hệ số này đ ợc tính từ tổng bình ph ơng trong kết qu phần (a) nh sau: R2 = 5.0139 + 2.1811 = 0.9074 5.0139 + 2.1811 + 0.7344 Và sau cùng, hệ số R2 điều chỉnh ph n ánh độ “c i tiến” của mô hình. Để hiểu hệ số này tốt hơn, chúng ta thấy ph ơng sai của toàn bộ mẫu là s2 = (5.0139 + 2.1811 + 0.7344) / 17 = 0.4644. Sau khi điều chỉnh cho nh h ng của điều kiện và vật liệu, ph ơng sai này còn 0.0525 (tức là residual mean square). Nh vậy hai yếu tố này làm gi m ph ơng sai kho ng 0.4644 – 0.0525 = 0.4119. Và hệ số R2 điều chỉnh là: Adj R2 = 0.4119 / 0.4644 = 0.88 Tức là sau khi điều chỉnh cho hai yếu tố điều kiện và vật liệu ph ơng sai của score gi m kho ng 88%. So sánh giữa các nhóm. Chúng ta sẽ ớc tính độ khác biệt giữa hai điều kiện và ba vật liệu bằng hàm TukeyHSD với aov: > res <- aov(score ~ condition+ material+condition) > TukeyHSD(res) Tukey multiple comparisons of means 95% family-wise confidence level 93 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Fit: aov(formula = score ~ condition + material + condition) $condition diff lwr upr p adj 2-1 -1.055556 -1.287131 -0.8239797 1e-07 $material diff lwr upr p adj 2-1 -0.8500000 -1.19610279 -0.5038972 0.0000442 3-1 -0.4833333 -0.82943612 -0.1372305 0.0068648 3-2 0.3666667 0.02056388 0.7127695 0.0374069 Biểu đồ sau đây sẽ minh ho cho các kết qu trên: > plot(TukeyHSD(res), ordered=TRUE) There were 16 warnings (use warnings() to see them) 3-2 3-1 2-1 95% family-wise confidence level -1.0 -0.5 0.0 0.5 Differences in mean levels of material Biểu đồ 25. So sánh giữa 3 lo i vật liệu bằng ph ơng pháp Tukey. 12. Phân tích hồi qui logistic Trong các phần tr ớc về phân tích hồi qui tuyến tính và phân tích ph ơng sai, chúng ta tìm mô hình và mối liên hệ giữa một biến phụ thuộc liên tục (continuous dependent variable) và một hay nhiều biến độc lập (independent variable) hoặc là liên tục hoặc là không liên tục. Nh ng trong nhiều tr ng hợp, biến phụ thuộc không ph i là biến liên tục mà là biến mang tính đo l ng nhị phân: có/không, mắc bệnh/không mắc bệnh, chết/sống, x y ra/không x y ra, v.v…, còn các biến độc lập có thể là liên tục hay không liên tục. Chúng ta cũng muốn tìm hiểu mối liên hệ giữa các biến độc lập và biến phụ thuộc. 94 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ví dụ 19. Trong một nghiên cứu do tôi tiến hành để tìm hiểu mối liên hệ giữa nguy cơ gãy x ơng (fracture, viết tắt là fx) và mật độ x ơng cùng một số chỉ số sinh hóa khác, 139 bệnh nhân nam (hay nói đúng hơn là đối t ợng nghiên cứu) tuổi từ 60 tr lên. Năm 1990, các số liệu sau đây đ ợc thu thập cho mỗi đối t ợng: độ tuổi (age), tỉ trọng cơ thể (body mass index hay BMI), mật độ chất khoáng trong x ơng (bone mineral density hay BMD), chỉ số hủy x ơng ICTP, chỉ số t o x ơng PINP. Các đối t ợng nghiên cứu đ ợc theo dõi trong vòng 15 năm. Trong th i gian theo dõi, các bệnh nhân bị gãy x ơng hay không gãy x ơng đ ợc ghi nhận. Câu hỏi đặt ra ban đầu là có một mối liên hệ gì giữa BMD và nguy cơ gãy x ơng hay không. Số liệu của nghiên cứu này đ ợc trình bày trong phần cuối của ch ơng này, và sẽ trình bày một phần d ới đây để b n đọc nắm đ ợc vấn đề. Một phần số liệu nghiên c u về các yếu tố nguy cơ cho gãy xương id 1 2 3 4 5 6 7 8 9 10 fx 1 1 1 1 1 0 0 0 0 0 age 79 89 70 88 85 68 70 69 74 79 bmi 24.7252 25.9909 25.3934 23.2254 24.6097 25.0762 19.8839 25.0593 25.6544 19.9594 bmd 0.818 0.871 1.358 0.714 0.748 0.935 1.040 1.002 0.987 0.863 ... ictp 9.170 7.561 5.347 7.354 6.760 4.939 4.321 4.212 5.605 5.204 pinp 37.383 24.685 40.620 56.782 58.358 67.123 26.399 47.515 26.132 60.267 137 138 139 0 1 0 64 80 67 38.0762 23.3887 25.9455 1.086 0.875 0.983 5.043 4.086 4.328 32.835 23.837 71.334 đây, vì biến phụ thuộc (gãy x ơng) không đ ợc đo l ng theo tính liên tục (mà chỉ là có hay không), cho nên ph ơng pháp phân tích hồi qui tuyến tính để phân tích mối liên hệ giữa biến phụ thuộc và biến độc lập. Một ph ơng pháp phân tích đ ợc phát triển t ơng đối gần đây (vào thập niên 1970s) có tên là logistic regression analysis (hay phân tích hồi qui logistic) có thể áp dụng cho tr ng hợp trên. Trong nghiên cứu này, sau 15 năm theo dõi, có 38 bệnh nhân bị gãy x ơng. Tính theo phần trăm, tỉ lệ gãy x ơng là 38 / 139 = 0.273 (hay 27.3%). 12.1 Mô hình hồi qui logistic Cho một tần số biến cố x ghi nhận từ n đối t ợng, chúng ta có thể tính xác suất của biến cố đó là: x p= n p có thể xem là một chỉ số đo l ng nguy cơ của một biến cố. Một cách thể hiện nguy cơ khác là odds (một danh từ, nếu tôi không lầm, chỉ có trong tiếng Anh – ngay c tiếng Pháp, Đức, Tây Ban Nha … cũng không có danh từ t ơng đ ơng với odds). Tôi t m dịch 95 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn odds là khả năng. Kh năng của một biến cố đ ợc định nghĩa đơn gi n bằng tỉ số xác suất biến cố x y ra trên xác suất biến cố không x y ra: p odds = 1− p Hàm logit của odds đ ợc định nghĩa nh sau:  p  l ogit ( p ) = log    1− p  Cho một biến độc lập x (x có thể là liên tục hay không liên tục), mô hình hồi qui logistic phát biểu rằng: logit(p) = α + β x T ơng tự nh mô hình hồi qui tuyến tính, α và β là hai thông số tuyến tính cần ph i ớc tính từ dữ liệu nghiên cứu. Nh ng ý nghĩa của thông số này, đặc biệt là thông số β, rất khác với ý nghĩa mà ta đã quen với mô hình hồi qui tuyến tính. Để hiểu ý nghĩa của hai thông số này, tôi sẽ quay l i với ví dụ 19. Vấn đề mà chúng ta muốn biết là mối liên hệ giữa mật độ x ơng bmd và nguy cơ gãy x ơng (fx). Để tiện cho việc minh họa, gọi bmd là x, vấn đề mà chúng ta cần biết có thể viết bằng ngôn ngữ mô hình nh sau  p  logit ( p ) = log  α + β x  1− p  Nói cách khác: odds ( p ) = p = eα + β x 1− p Nói cách khác, mô hình hồi qui logistic vừa trình bày trên phát biểu rằng mối liên hệ giữa xác suất gãy x ơng (p) và mật độ x ơng bmd là một mối liên hệ theo hình chữ S. Mô hình trên còn cho thấy xác suất gãy x ơng p tùy thuộc vào giá trị của x. Thành ra, mô hình trên có thể viết một cách chính xác hơn rằng khả năng gãy x ơng với điều kiện x là: odds ( p | x ) = eα + β x Khi x = x0, kh năng gãy x ơng là: odds ( p | x = x0 ) = eα + β x0 Khi x = x0 + 1 (tức tăng 1 đơn vị từ x0), kh năng gãy x ơng là: odds ( p | x = x0 + 1) = e α + β ( x0 +1) Và, tỉ số của hai xác suất gãy x ơng: 96 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn odds ( p | x = x0 + 1) odds ( p | x = x0 ) = α + β ( x0 +1) e = eβ α + β x0 e Trong dịch tễ học, eβ đ ợc gọi là odds ratio. Odds ratio, nh tên gọi là, tỉ số khả năng hay tỉ số khả dĩ. Nói cách khác, hệ số β trong mô hình hồi qui logistic chính là tỉ số kh dĩ. Ph ơng pháp để ớc tính thông số trong mô hình [3] khá phức t p (dùng ph ơng pháp maximum likelihood – tức ph ơng pháp Hợp lí cực đại) và không nằm trong ph m vi của cuốn sách này, nên tôi sẽ không trình bày đây (b n đọc có thể tham kh o sách giáo khoa để biết thêm, nếu cần thiết). Tuy nhiên, tôi muốn đề cập ngắn gọn là ph ơng pháp hợp lí cực đ i cung cấp cho chúng ta một hệ ph ơng trình nh sau: ( ) −1 n  n −(αˆ + βˆ xi ) ∑ yi = ∑ 1 + e  i =1 i =1  n n  x y = x 1 + e −(αˆ + βˆ xi ) ∑ i i i ∑ i =1 i =1 ( ) Trong đó, Trong đó, yi là biến phụ thuộc (gãy x ơng với giá trị 0 hay 1), và xi là biến độc lập (mật độ x ơng), và n là số mẫu. Để tìm ớc số α̂ và βˆ , một trong những phép tính hay sử dụng là iterative weighted least square hay Newton-Raphson. R sử dụng phép tính Newton-Raphson để tìm hai ớc số đó. Sau khi đã có ớc số α̂ và βˆ chúng ta có thể ớc tính xác suất p cho bất cứ giá trị nào của x nh sau (sau vài thao tác đ i số): pˆ = eαˆ + β x 1+ e ˆ αˆ + βˆ x = 1+ e 1 ( − αˆ + βˆ x ) Chú ý tôi dùng dấu mũ p̂ để chỉ số ớc tính (predicted value), chứ không ph i p là xác suất quan sát. Nếu mô hình mô t dữ liệu tốt và đầy đủ, độ khác biệt giữa p và p̂ nhỏ; nếu mô hình không thích hợp hay không tốt, độ khác biệt đó có thể sẽ cao. Độ khác biệt giữa p và p̂ đ ợc gọi là deviance. Ph ơng pháp tính deviance khá phức t p, nh ng đó không ph i là chủ đề đây, cho nên tôi chỉ nói qua khái niệm mà thôi. Khi chúng ta có nhiều mô hình để mô phỏng một hay nhiều mối liên hệ, deviance có thể đ ợc sử dụng để đánh giá sự thích hợp của một mô hình, hay để chọn một mô hình “tối u”. 12.2 Phân tích hồi qui logistic bằng R Bây gi , chúng ta quay l i với ví dụ 1, dùng số liệu trong B ng 12.1 để ớc tính hai thông số α và β bằng R. Tr ớc hết chúng ta ph i nhập toàn bộ số liệu vào một data 97 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn frame, và cho một cái tên, chẳng h n nh fracture. Trong tr ng hợp của tôi, dữ liệu đ ợc chứa trong directory c:\works\stats d ới tên fracture.txt, do đó, các lệnh sau đây cần thiết để nhập số liệu: # báo cho R biết nơi chứa số liệu > setwd(“c:/works/stats”) # nhập số liệu và cho vào một data frame tên fracture > fracture <- read.table(“fracture.txt”, header=TRUE, na.string=”.”) # kiểm tra xem có bao nhiêu biến trong dữ liệu fracture > names(fracture) [1] "id" "fx" "age" "bmi" "bmd" "ictp" "pinp" # Chọn những bệnh nhân có đầy đủ số liệu cho phân tích > fulldata <- na.omit(fracture) > attach(fulldata) Hai biến mà chúng ta quan tâm trong ví dụ này là: fx (gãy x ơng) và bmd (mật độ x ơng). Chúng ta kiểm tra xem có bao nhiêu bệnh nhân gãy x ơng: > table(fx) fx 0 1 101 38 Kế đến, xem mật độ x ơng trong nhóm gãy x ơng và không gãy x ơng ra sao: > tapply(bmd, fx, mean) 0 1 0.9444851 0.9016667 > boxplot(bmd ~ fx, xlab=”Fracture: 1=yes, 0=no)”, ylab=”BMD”) 98 Phân tích số liệu và biểu đồ bằng R 1.0 0.6 0.8 BMD 1.2 Nguyễn Văn Tuấn 0 1 Fracture: 1=yes, 0=no) Kết qu trên cho thấy, bmd trong nhóm bệnh nhân bị gãy x ơng thấp hơn so với nhóm không bị gãy x ơng (0.90 và 0.94). Và, kiểm định t sau đây cho thấy mức độ khác biệt này không có ý nghĩa thống kê (p = 0.15). > t.test(bmd~fx) Welch Two Sample t-test data: bmd by fx t = 1.4572, df = 53.952, p-value = 0.1508 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.01609226 0.10172922 sample estimates: mean in group 0 mean in group 1 0.9444851 0.9016667 Để ớc tính thông số trong mô hình [4], hàm số glm (viết tắt từ generalized linear model) trong R có thể áp dụng, với “cú pháp” nh sau: > logistic <- glm(fx ~ bmd, family=”binomial”) > summary(logistic) Call: glm(formula = fx ~ bmd, family = "binomial") Deviance Residuals: Min 1Q Median -1.0287 -0.8242 -0.7020 3Q 1.3780 Max 2.0709 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.063 1.342 0.792 0.428 bmd -2.270 1.455 -1.560 0.119 (Dispersion parameter for binomial family taken to be 1) 99 Phân tích số liệu và biểu đồ bằng R Null deviance: 157.81 Residual deviance: 155.27 AIC: 159.27 Nguyễn Văn Tuấn on 136 on 135 degrees of freedom degrees of freedom Number of Fisher Scoring iterations: 4 Tôi sẽ lần l ợt gi i thích các kết qu trên: (a) Trong lệnh logistic <- glm(fx ~ bmd, family=”binomial”) chúng ta yêu cầu R phân tích theo mô hình fx là một hàm số với bmd nh mô hình [4]. Trong glm có nhiều luật phân phối, mà trong đó phân phối nhị phân (binomial) là một luật phân phối chuẩn cho hồi qui logistic. Do đó, family=”binomial” cần thiết cho R. (b) Deviance: phần thứ nhất của kết qu cho biết qua về deviance. Deviance Residuals: Min 1Q Median -1.0287 -0.8242 -0.7020 3Q 1.3780 Max 2.0709 Deviance nh gi i thích trên ph n ánh độ khác biệt giữa mô hình và dữ liệu (cũng t ơng tự nh mean square residual trong phân tích hồi qui tuyến tính vậy). Đối với một mô hình đơn lẻ nh ví dụ này thì giá trị của deviance không có ý nghĩa gì nhiều. (c) Phần kế tiếp cung cấp ớc số của α̂ (mà R đặt tên là intercept) và βˆ (bmd) và sai số chuẩn (standard error). Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.063 1.342 0.792 0.428 bmd -2.270 1.455 -1.560 0.119 Qua kết qu này, chúng ta có α̂ = 1.063 và βˆ = -2.27. ớc số βˆ là số âm cho thấy mối liên hệ giữa nguy cơ gãy x ơng và bmd là mối liên hệ nghịch đ o: xác suất gãy x ơng tăng khi giá trị của bmd gi m. Tuy nhiên, kiểm định z (tính bằng cách lấy ớc số chia cho sai số chuẩn) cho chúng ta thấy nh h ng của bmd không có ý nghĩa thống kê, vì trị số p = 0.119. Nhớ rằng tỉ số kh dĩ (odds ratio hay viết tắt là OR) chính là e-2.27 = 0.1033. Nói cách khác, khi bmd tăng 1 g/cm2 (đơn vị đo l ng của bmd là g/cm2) thì tỉ số OR gi m 0.9067 hay 90.67%. Nh ng tăng 1 g/cm2 là mật độ rất cao trong x ơng và không thực tế. Cho nên một cách tính khác là tính trên độ lệch chuẩn (standard deviation) của bmd. Chúng ta sẽ tìm hiểu độ lệch chuẩn của bmd: > sd(bmd) [1] 0.1406543 Do đó, OR sẽ tính trên mỗi 0.14 g/cm2. Và OR cho mỗi độ lệch chuẩn, do đó, là: 100 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn e-2.27*0.1406 = 0.7267 Tức là, khi bmd tăng một độ lệch chuẩn thì tỉ số kh dĩ gãy x ơng gi m kho ng 28%. Cũng có thể nói cách khác, là khi bmd giảm một độ lệch chuẩn thì tỉ số kh dĩ tăng e2.27*0.1406 = 1.376 hay kho ng 38%. Một cách khác để biết nh h trình: ng của bmd là ớc tính xác suất gãy x ơng qua ph ơng pˆ = ( ) e 1.063− 2.27 ( bmd ) 1+ e 1.063− 2.27 bmd Theo đó, khi bmd = 1.00, p = 0.23. Khi bmd = 0.86 (tức gi m 1 độ lệch chuẩn), p = 0.291. Tức là, nếu BMD gi m 1 độ lệch chuẩn thì xác suất gãy x ơng tăng 0.291/0.23 = 1.265 hay 26%5. (d) Phần cuối của kết qu cung cấp deviance cho hai mô hình: mô hình không có biến độc lập (null deviance), và mô hình với biến độc lập, tức là bmd trong ví dụ (residual deviance). Null deviance: 157.81 Residual deviance: 155.27 AIC: 159.27 on 136 on 135 degrees of freedom degrees of freedom Qua hai số này, chúng ta thấy bmd nh h ng rất thấp đến việc tiên đoán gãy x ơng, chỉ làm gi m deviance từ 157.8 xuống còn 155.27, và mức độ gi m này không có ý nghĩa thống kê. Ngoài ra, R còn cung cấp giá trị của AIC (Akaike Information Criterion) đ ợc tính từ deviance và bậc tự do. Tôi sẽ quay l i ý nghĩa của AIC trong phần sắp đến khi so sánh các mô hình. 12.3 ớc tính xác suất bằng R Xin nhắc l i trong phân tích trên, chúng ta cho các kết qu vào đối t ợng logistic. Trong đối t ợng này có nhiều thông tin có ích, nh ng nếu muốn xem các thông tin này chúng ta ph i dùng đến các lệnh nh summary chẳng h n. Trong phần này, tôi sẽ trình bày một vài hàm để xem xét các thông tin liên quan đến việc tiên đoán xác suất. • predict dùng để liệt kê các giá trị ớc tính (predicted values) của mô hình  p  log   = α + β x cho từng bệnh nhân.  1− p  > predict(logistic) 1 2 3 101 4 5 6 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 2.377576584 1.085694014 -2.141117756 1.492824115 0.965379946 -0.941253280 7 8 9 10 11 12 -1.733686514 -1.675645430 -0.665282957 -0.507046129 -0.941854868 -0.648740461 ... Các số trên là log(p / (1 – p)), tức log odds, không có ý nghĩa hực tế bao nhiêu. Chúng ta 1.063− 2.27 ( bmd ) e muốn biết giá trị tiên đoán xác suất p tính từ ph ơng trình pˆ = . Để có giá trị 1.063− 2.27 ( bmd ) 1+ e này cho từng bệnh nhân, chúng ta cho thông số type=”response” vào hàm predict nh sau: > predict(logistic, type="response") 1 2 3 4 5 6 7 0.91510135 0.74757001 0.10516416 0.81650178 0.72419767 0.28064726 0.15011664 8 9 10 11 12 13 14 0.15767295 0.33955387 0.37588624 0.28052582 0.34327343 0.44305196 0.23830776 ... Trong kết qu trên (chỉ in một phần) ớc tính xác suất gãy x ơng cho bệnh nhân 1 là 0.915, cho bệnh nhân 2 là 0.747, v.v… • Chúng ta có thể xem xét các giá trị tiên đoán này với độ bmd bằng cách dùng hàm plot thông th ng: 0.35 0.30 0.25 0.20 0.15 fitted(glm(fx ~ bmd, family = "binomial")) 0.40 > plot(bmd, fitted(glm(fx ~ bmd, family=”binomial”))) 0.6 0.8 1.0 1.2 bmd Xác suất tiên đoán gãy x ơng (trục tung) và độ bmd (trục hoành) qua mô hình hồi qui logistic. 102 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Biểu đồ trên có thể c i tiến bằng cách cho các kho ng cách giá trị bmd gần nhau hơn (nh 0.50, 0.55, 0.60, …, 1.20 chẳng h n), và dùng đ ng biểu diễn thay vì dùng dấu chấm. Các lệnh sau đây sẽ c i tiến biểu đồ. logistic <- glm(fx ~ bmd, family=”binomial”) fnbmd <- seq(0.5, 1.2, 0.05) #cho fnbmd từ > 0.50,0.55,0.6,...,1.2 new.data <- data.frame(bmd = fnbmd) #cho vào một dataframe mới predicted <- predict(logistic, new.data, type=”response”) plot(predicted ~ fnbmd, type=”l”) 0.35 0.30 0.15 0.20 0.25 predicted 0.40 0.45 > > > > > 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 fnbmd Xác suất tiên đoán gãy x ơng (trục tung) và độ bmd (trục hoành) qua mô hình hồi qui logistic. 13. ớc tính cỡ mẫu (sample size estimation) Một công trình nghiên cứu th ng dựa vào một mẫu (sample). Một trong những câu hỏi quan trọng nhất tr ớc khi tiến hành nghiên cứu là cần bao nhiêu mẫu hay bao nhiêu đối t ợng cho nghiên cứu. “Đối t ợng” đây là đơn vị căn b n của một nghiên cứu, là số bệnh nhân, số tình nguyện viên, số mẫu ruộng, cây trồng, thiết bị, v.v… ớc tính số l ợng đối t ợng cần thiết cho một công trình nghiên cứu đóng vai trò cực kì quan trọng, vì nó có thể là yếu tố quyết định sự thành công hay thất b i của nghiên cứu. Nếu số l ợng đối t ợng không đủ thì kết luận rút ra từ công trình nghiên cứu không có độ chính xác cao, thậm chí không thể kết luận gì đ ợc. Ng ợc l i, nếu số l ợng đối t ợng quá nhiều hơn số cần thiết thì tài nguyên, tiền b c và th i gian sẽ bị hao phí. Do đó, vấn đề then chốt tr ớc khi nghiên cứu là ph i ớc tính cho đ ợc một số đối t ợng vừa đủ cho mục tiêu của nghiên cứu. Số l ợng đối t ợng “vừa đủ” tùy thuộc vào ba yếu tố chính: 103 Phân tích số liệu và biểu đồ bằng R • • • Nguyễn Văn Tuấn Sai sót mà nhà nghiên cứu chấp nhận, cụ thể là sai sót lo i I và II; Độ dao động (variability) của đo l ng, mà cụ thể là độ lệch chuẩn; và Mức độ khác biệt hay nh h ng mà nhà nghiên cứu muốn phát hiện. Không có số liệu về ba yếu tố này thì không thể nào ớc tính cỡ mẫu. Kinh nghiệm của ng i viết cho thấy rất nhiều ng i khi tiến hành nghiên cứu th ng không có ý niệm gì về các số liệu này, cho nên khi đến tham vấn các chuyên gia về thống kê học, họ chỉ nhận câu tr l i: “không thể tính đ ợc”! Trong ch ơng này tôi sẽ bàn qua ba yếu tố trên. 13.1 Khái niệm về “power” Thống kê học là một ph ơng pháp khoa học có mục đích phát hiện, hay đi tìm những cái có thể gộp chung l i bằng cụm từ “ch a đ ợc biết” (unknown). Cái ch a đ ợc biết đây là những hiện t ợng chúng ta không quan sát đ ợc, hay quan sát đ ợc nh ng không đầy đủ. “Cái ch a biết” có thể là một ẩn số (nh chiều cao trung bình ng i Việt Nam, hay trọng l ợng một phần tử), hiệu qu của một thuật điều trị, gen có chức năng làm cho cây lá có màu xanh, s thích của con ng i, v.v… Chúng ta có thể đo chiều cao, hay tiến hành xét nghiệm để biết hiệu qu của thuốc, nh ng các nghiên cứu nh thế chỉ đ ợc tiến hành trên một nhóm đối t ợng, chứ không ph i toàn bộ quần thể của dân số. mức độ đơn gi n nhất, những cái ch a biết này có thể xuất hiện d ới hai hình thức: hoặc là có, hoặc là không. Chẳng h n nh một thuật điều trị có hay không có hiệu qu chống gãy x ơng, khách hàng thích hay không thích một lo i n ớc gi i khát. B i vì không ai biết hiện t ợng một cách đầy đủ, chúng ta ph i đặt ra gi thiết. Gi thiết đơn gi n nhất là giả thiết đảo (hiện t ợng không tồn t i, kí hiệu H-) và giả thiết chính (hiện t ợng tồn t i, kí hiệu H+). Chúng ta sử dụng các ph ơng pháp kiểm định thống kê (statistical test) nh kiểm định t, F, z, χ2, v.v… để đánh giá kh năng của gi thiết. Kết qu của một kiểm định thống kê có thể đơn gi n chia thành hai giá trị: hoặc là có ý nghĩa thống kê (statistical significance), hoặc là không có ý nghĩa thống kê (non-significance). Có ý nghĩa thống kê đây, nh đề cập trong Ch ơng 7, th ng dựa vào trị số P: nếu P < 0.05, chúng ta phát biểu kết qu có ý nghĩa thống kê; nếu P > 0.05 chúng ta nói kết qu không có ý nghĩa thống kê. Cũng có thể xem có ý nghĩa thống kê hay không có ý nghĩa thống kê nh là có tín hiệu hay không có tín hiệu. Hãy t m đặt kí hiệu T+ là kết qu có ý nghĩa thống kê, và T- là kết qu kiểm định không có ý nghĩa thống kê. Hãy xem xét một ví dụ cụ thể: để biết thuốc risedronate có hiệu qu hay không trong việc điều trị loãng x ơng, chúng ta tiến hành một nghiên cứu gồm 2 nhóm bệnh nhân (một nhóm đ ợc điều trị bằng risedronate và một nhóm chỉ sử dụng gi d ợc placebo). Chúng ta theo dõi và thu thập số liệu gãy x ơng, ớc tính tỉ lệ gãy x ơng cho từng nhóm, và so sánh hai tỉ lệ bằng một kiểm định thống kê. Kết qu kiểm định thống kê hoặc là có ý nghĩa thống kê (P<0.05) hay không có ý nghĩa thống kê (P>0.05). Xin nhắc l i rằng chúng ta không biết risedronate thật sự có hiệu nghiệm chống gãy x ơng 104 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn hay không; chúng ta chỉ có thể đặt gi thiết H. Do đó, khi xem xét một gi thiết và kết qu kiểm định thống kê, chúng ta có bốn tình huống: (a) Gi thuyết H đúng (thuốc risedronate có hiệu nghiệm) và kết qu kiểm định thống kê P<0.05. (b) Gi thuyết H đúng, nh ng kết qu kiểm định thống kê không có ý nghĩa thống kê; (c) Gi thuyết H sai (thuốc risedronate không có hiệu nghiệm) nh ng kết qu kiểm định thống kê có ý nghĩa thống kê; (d) Gi thuyết H sai và kết qu kiểm định thống kê không có ý nghĩa thống kê. đây, tr ng hợp (a) và (d) không có vấn đề, vì kết qu kiểm định thống kê nhất quán với thực tế của hiện t ợng. Nh ng trong tr ng hợp (b) và (c), chúng ta ph m sai lầm, vì kết qu kiểm định thống kê không phù hợp với gi thiết. Trong ngôn ngữ thống kê học, chúng ta có vài thuật ngữ: • • • • xác suất của tình huống (b) x y ra đ ợc gọi là sai sót loại II (type II error), và th ng kí hiệu bằng β. xác suất của tình huống (a) đ ợc gọi là Power. Nói cách khác, power chính là xác suất mà kết qu kiểm định thống cho ra kết qu p<0.05 với điều kiện gi thiết H là thật. Nói cách khác: power = 1-β ; xác suất của tình huống (c) đ ợc gọi là sai sót loại I (type I error, hay significance level), và th ng kí hiệu bằng α. Nói cách khác, α chính là xác suất mà kết qu kiểm định thống cho ra kết qu p<0.05 với điều kiện gi thiết H sai; xác suất tình hống (d) không ph i là vấn đề cần quan tâm, nên không có thuật ngữ, dù có thể gọi đó là kết qu âm tính thật (hay true negative). Có thể tóm l ợc 4 tình huống đó trong một B ng 1 sau đây: Các tình huống trong việc thử nghiệm một giả thiết khoa học Kết quả kiểm định thống kê Đúng (thuốc có hiệu nghiệm) Giả thuyết H Sai (thuốc không có hiệu nghiệm) Có ý nghĩa thống kê (p<0,05) D ơng tính thật (power), 1-β= P(s | H+) Sai sót lo i I (type I error) α = P(s | H-) Không có ý nghĩa thống kê (p>0,05) Sai sót lo i II (type II error) β = P(ns | H+) Âm tính thật (true negative) 1-α = P(ns | H-) 105 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Chú thích: s trong biểu đồ này có nghĩa là significant; ns non-significant; H+ là gi thuyết đúng; và H- là gi thuyết sai. Do đó, có thể mô t 4 tình huống trên bằng ngôn ngữ xác suất có điều kiện nh sau: Power = 1 – β = P(s | H+); β = P(ns | H+); và α = P(s | H-). 13.2 Số liệu để ớc tính cỡ mẫu Nh đã đề cập trong phần đầu của ch ơng này, để ớc tính số đối t ợng cần thiết cho một công trình nghiên cứu, chúng ta cần ph i có 3 số liệu: xác suất sai sót lo i I và II, độ dao động của đo l ng, và độ nh h ng. • • • Về xác suất sai sót, thông th ng một nghiên cứu chấp nhận sai sót lo i I kho ng 1% hay 5% (tức α = 0.01 hay 0.05), và xác suất sai sót lo i II kho ng β = 0.1 đến β = 0.2 (tức power ph i từ 0.8 đến 0.9). Độ dao động chính là độc lệch chuẩn (standard deviation) của đo l ng mà công trình nghiên cứu dựa vào để phân tích. Chẳng h n nh nếu nghiên cứu về cao huyết áp, thì nhà nghiên cứu cần ph i có độ lệch chuẩn của áp suất máu. Chúng ta t m gọi độ dao động là σ. Độ nh h ng, nếu là công trình nghiên cứu so sánh hai nhóm, là độ khác biệt trung bình giữa hai nhóm mà nhà nghiên cứu muốn phát hiện. Chẳng h n nh nhà nghiên cứu có thể gi thiết rằng bệnh nhân đ ợc điều trị bằng thuốc A có àp suất máu gi m 10 mmHg so với nhóm gi đ ợc. đây, 10 mmHg đ ợc xem là độ nh h ng. Chúng ta t m gọi độ nh h ng là ∆. Một nghiên cứu có thể có một nhóm đối t ợng hay hai (và có khi hơn 2) nhóm đối t ợng. Và ớc tính cỡ mẫu cũng tùy thuộc vào các tr ng hợp này. Trong tr ng hợp một nhóm đối t ợng, số l ợng đối t ợng (n) cần thiết cho nghiên cứu có thể tính toán một cách “thủ công” nh sau: n= (∆ /σ ) C 2 Trong tr ng hợp có hai nhóm đối t ợng, số l ợng đối t ợng (n) cần thiết cho nghiên cứu có thể tính toán nh sau: n = 2× (∆ /σ ) C 2 Trong đó, hằng số C đ ợc xác định từ xác suất sai sót lo i I và II (hay power) nh sau: 106 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Hằng số C liên quan đến sai sót loại I và II α= β = 0.20 (Power = 0.80) 6.15 7.85 13.33 0.10 0.05 0.01 β = 0.10 (Power = 0.90) 8.53 10.51 16.74 β = 0.05 (Power = 0.95) 10.79 13.00 19.84 ớc tính cỡ mẫu 13.4 13.4.1 ớc tính cỡ mẫu cho một chỉ số trung bình Ví dụ 20: Chúng ta muốn ớc tính chiều cao đàn ông ng i Việt, và chấp nhận sai số trong vòng 1 cm (d = 1) với kho ng tin cậy 0.95 (tức α=0.05) và power = 0.8 (hay β = 0.2). Các nghiên cứu tr ớc cho biết độ lệch chuẩn chiều cao ng i Việt kho ng 4.6 cm. Chúng ta có thể áp dụng công thức [1] để ớc tính cỡ mẫu cần thiết cho nghiên cứu: n= (∆ /σ ) C 2 = (1/ 4.6 ) Nói cách khác, chúng ta cần ph i đo chiều cao ông Việt với sai số trong vòng 1 cm. 7.85 2 = 166 166 đối t ợng để ớc tính chiều cao đàn Nếu sai số chấp nhận là 0.5 cm (thay vì 1 cm), số l ợng đối t ợng cần thiết là: 7.85 n= = 664 . Nếu độ sai số mà chúng ta chấp nhận là 0.1 cm thì số l ợng đối 2 ( 0.5 / 4.6 ) t ợng nghiên cứu lên đến 16610 ng i! Qua các ớc tính này, chúng ta dễ dàng thấy cỡ mẫu tùy thuộc rất lớn vào độ sai số mà chúng ta chấp nhận. Muốn có ớc tính càng chính xác, chúng ta cần càng nhiều đối t ợng nghiên cứu. Trong R có hàm power.t.test có thể áp dụng để ớc tính cỡ mẫu cho ví dụ trên nh sau. Chú ý chúng ta cho R biết vấn đề là một nhóm tức type=”one.sample”: # sai số 1 cm, độc lệch chuẩn 4.6, a=0.05, power=0.8 > power.t.test(delta=1, sd=4.6, sig.level=.05, power=.80, type='one.sample') One-sample t test power calculation n delta sd sig.level power alternative = = = = = = 168.0131 1 4.6 0.05 0.8 two.sided 107 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn kết qu tính toán từ R là 168, khác với cách tính thủ công 2 đối t ợng, vì cố nhiên R sử dụng nhiều số lẻ hơn và chính xác hơn cách tính thủ công. Với sai số 0.5 cm: # sai số 0.5 cm, độc lệch chuẩn 4.6, a=0.05, power=0.8 > power.t.test(delta=0.5, sd=4.6, sig.level=.05, power=.80, type='one.sample') One-sample t test power calculation n delta sd sig.level power alternative = = = = = = 666.2525 0.5 4.6 0.05 0.8 two.sided Ví dụ 21: Một lo i thuốc điều trị có kh năng tăng độ alkaline phosphatase bệnh nhân loãng x ơng. Độ lệch chuẩn của alkaline phosphatase là 15 U/l. Một nghiên cứu mới sẽ tiến hành trong một quần thể bệnh nhân Việt Nam, và các nhà nghiên cứu muốn biết bao nhiêu bệnh nhân cần tuyển để chứng minh rằng thuốc có thể alkaline phosphatase từ 60 đến 65 U/l sau 3 tháng điều trị, với sai số I α = 0.05 và power = 0.8. Đây là một lo i nghiên cứu “tr ớc – sau” (before-after study); có nghĩa là tr ớc và sau khi điều trị. đây, chúng ta chỉ có một nhóm bệnh nhân, nh ng đ ợc đo hai lần (tr ớc khi dùng thuốc và sau khi dùng thuốc). Chỉ tiêu lâm sàng để đánh giá hiệu nghiệm của thuốc là độ thay đổi về alkaline phosphatase. Trong tr ng hợp này, chúng ta có trị số tăng trung bình là 5 U/l và độ lệch chuẩn là 15 U/l, hay nói theo ngôn ngữ R, delta=5, sd=15, sig.level=.05, power=.80, và lệnh: > power.t.test(delta=3, sd=15, sig.level=.05, power=.80, type='one.sample') One-sample t test power calculation n delta sd sig.level power alternative = = = = = = 198.1513 3 15 0.05 0.8 two.sided Nh vậy, chúng ta cần ph i có 198 bệnh nhân để đ t các mục tiêu trên. 13.4.2 ớc tính cỡ mẫu cho so sánh hai số trung bình Trong thực tế, rất nhiều nghiên cứu nhằm so sánh hai nhóm với nhau. Cách ớc tính cỡ mẫu cho các nghiên cứu này chủ yếu dựa vào công thức [2] nh trình bày phần 15.3.1. 108 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ví dụ 22: Một nghiên cứu đ ợc thiết kế để thử nghiệm thuốc alendronate trong việc điều trị loãng x ơng phụ nữ sau th i kì mãn kinh. Có hai nhóm bệnh nhân đ ợc tuyền: nhóm 1 là nhóm can thiệp (đ ợc điều trị bằng alendronate), và nhóm 2 là nhóm đối chứng (tức không đ ợc điều trị). Tiêu chí để đánh giá hiệu qu của thuốc là mật độ x ơng (bone mineral density – BMD). Số liệu từ nghiên cứu dịch tễ học cho thấy giá trị trung bình của BMD trong phụ nữ sau th i kì mãn kinh là 0.80 g/cm2, với độ lệch chuẩn là 0.12 g/cm2. Vấn đề đặt ra là chúng ta cần ph i nghiên cứu bao nhiêu đối t ợng để “chứng minh” rằng sau 12 tháng điều trị BMD của nhóm 1 tăng kho ng 5% so với nhóm 2? Trong ví dụ trên, t m gọi trị số trung bình của nhóm 2 là µ 2 và nhóm 1 là µ1 , chúng ta có: µ1 = 0.8*1.05 = 0.84 g/cm2 (tức tăng 5% so với nhóm 1), và do đó, ∆ = 0.84 – 0.80 = 0.04 g/cm2. Độ lệch chuẩn là σ = 0.12 g/cm2. Với power = 0.90 và α = 0.05, cỡ mẫu cần thiết là: n= (∆ /σ ) 2C 2 = 2 ×10.51 ( 0.04 / 0.12 ) 2 = 189 Và l i gi i từ R qua hàm power.t.test nh sau: > power.t.test(delta=0.04, sd=0.12, sig.level=0.05, power=0.90, type="two.sample") Two-sample t test power calculation n delta sd sig.level power alternative = = = = = = 190.0991 0.04 0.12 0.05 0.9 two.sided NOTE: n is number in *each* group Chú ý trong hàm power.t.test, ngoài các thông số thông th ng nh delta (độ nh h ng hay khác biệt theo gi thiết), sd (độ lệch chuẩn), sig.level xác suất sai sót lo i I, và power, chúng ta còn ph i cụ thể chỉ ra rằng đây là nghiên cứu gồm có hai nhóm với thông số type=”two.sample”. Kết qu trên cho biết chúng ta cần 190 bệnh nhân cho mỗi nhóm (hay 380 bệnh nhân cho công trình nghiên cứu). Trong tr ng hợp này, power = 0.90 và α = 0.05 có nghĩa là gì ? Tr l i: hai thông số đó có nghĩa là nếu chúng ta tiến hành thật nhiều nghiên cứu (ví dụ 1000) và mỗi nghiên cứu với 380 bệnh nhân, sẽ có 90% (hay 900) nghiên cứu sẽ cho ra kết qu trên với trị số p < 0.05. 109 Phân tích số liệu và biểu đồ bằng R 13.4.3 Nguyễn Văn Tuấn ớc tính cỡ mẫu cho phân tích ph ơng sai Ph ơng pháp ớc tính cỡ mẫu cho so sánh giữa hai nhóm cũng có thể khai triển thêm để ớc tính cỡ mẫu cho tr ng hợp so sánh hơn hai nhóm. Trong tr ng hợp có nhiều nhóm, nh đề cập trong Ch ơng 11, ph ơng pháp so sánh là phân tích ph ơng sai. Theo ph ơng pháp này, số trung bình bình ph ơng phần d (residual mean square, RMS) chính là ớc tính của độ dao động của đo l ng trong mỗi nhóm, và chỉ số này rất quan trọng trong việc ớc tính cỡ mẫu. Chi tiết về lí thuyết đằng sau cách ớc tính cỡ mẫu cho phân tích ph ơng sai khá phức t p, và không nằm trong ph m vi của ch ơng này. Nh ng nguyên lí chủ yếu vẫn không khác so với lí thuyết so sánh giữa hai nhóm. Gọi số trung bình của k nhóm là µ1, µ2, µ3, . . ., µk, chúng ta có thể tính tổng bình ph ơng giữa các nhóm bằng k k SS 2 SS SS = ∑ ( µi − µ ) , trong đó, µ = ∑ µi / k . Cho λ = , vấn đề đặt ra là tìm ( k − 1) RMS i =1 i =1 cố l ợng cỡ mẫu n sao cho zβ đáp ứng yêu cầu power = 0.80 hay 0.9, mà zβ = ( k − 1)(1 + nλ ) F + k ( n − 1)(1 + 2nλ ) 1 × 2    k ( n − 1)  2 ( k − 1)(1 + nλ ) − (1| 2nλ )  − F ( k − 1)(1 + nλ ) ( 2k ( n − 1) − 1)    Trong đó F là kiểm định F. (Xem J. Fleiss, “The Design and Analysis of Clinical Experiments”, John Wiley & Sons, New York 1986, trang 373). Ví dụ 23. Để so sánh độ ngọt của một lo i n ớc uống giữa 4 nhóm đối t ợng khác nhau về giới tính và độ tuổi (t m gọi 4 nhóm là A, B, C và D), các nhà nghiên cứu gi thiết rằng độ ngọt trong nhóm A, B. C và D lần l ợc là 4.5, 3.0, 5.6, và 1.3. Qua xem xét nhiều nghiên cứu tr ớc, các nhà nghiên cứu còn biết rằng RMS về độ ngọt trong mỗi nhóm là kho ng 8.7. Vấn đề đặt ra là bao nhiêu đối t ợng cần nghiên cứu để phát hiện sự khác biệt có ý nghĩa thống kê mức độ α = 0.05 và power = 0.9. Hàm power.anova.test trong R có thể ứng dụng để gi i quyết vấn đề. Chúng ta chỉ cần đơn gi n cung cấp 4 số trung bình theo gi thiết và số RMS nh sau: # trước hết cho 4 số trung bình vào một vector > groupmeans <- c(4.5, 3.0, 5.6, 1.3) # sau đó, “gọi” hàm power.anova.test: > power.anova.test(groups = length(groupmeans), between.var=var(groupmeans), within.var=8.7, power=0.90, sig.level=0.05) Balanced one-way analysis of variance power calculation groups = 4 110 Phân tích số liệu và biểu đồ bằng R n between.var within.var sig.level power = = = = = Nguyễn Văn Tuấn 12.81152 3.486667 8.7 0.05 0.9 NOTE: n is number in each group Kết qu cho thấy các nhà nghiên cứu cần kho ng 13 đối t ợng cho mỗi nhóm (tức 52 đối t ợng cho toàn bộ nghiên cứu). 13.4.4 ớc tính cỡ mẫu để ớc tính một tỉ lệ Nhiều nghiên cứu mô t có mục đích khá đơn gi n là ớc tính một tỉ lệ. Chẳng h n nh giới y tế th ng hay tìm hiểu tỉ lệ một bệnh trong cộng đồng, hay giới thăm dò ý kiến và thị tr ng th ng tìm hiểu tỉ lệ dân số a thích một s n phẩm. Trong các tr ng hợp này, chúng ta không có những đo l ng mang tính liên tục, nh ng kết qu chỉ là những giá trị nhị nh có / không, thích / không tích, v.v… Và cách ớc tính cỡ mẫu cũng khác với ba ví dụ trên đây. Năm 1991, một cuộc thăm dò ý kiến Mĩ cho thấy 45% ng i đ ợc hỏi sẵn sàng khuyến khích con họ nên hiến một qu thận cho những bệnh nhân cần thiết. Kho ng tin cậy 95% của tỉ lệ này là 42% đến 48%, tức một kho ng cách đến 6%! Kết qu này [t ơng đối] thiếu chính xác, dù số l ợng đối t ợng tham gia lên đến 1000 ng i. T i sao? Để tr l i câu hỏi này, chúng ta thử xem qua một vài lí thuyết về ớc tính cỡ mẫu cho một tỉ lệ. Chúng ta biết qua Ch ơng 6 và 9 rằng nếu p̂ đ ợc ớc tính từ n đối t ợng, thì kho ng tin cậy 95% của một tỉ lệ p [trong dân số] là: pˆ ± 1.96 × SE ( pˆ ) , trong đó SE ( pˆ ) = pˆ (1 − pˆ ) / n . Bây gi thử lật ng ợc vấn đề: chúng ta muốn ớc tính p sao kho ng rộng 2 × 1.96 × SE ( pˆ ) không quá một hằng số m. Nói cách khác, chúng ta muốn: 1.96 × pˆ (1 − pˆ ) / n ≤ m Chúng ta muốn tìm số l ợng đối t ợng n để đ t yêu câu trên. Qua cách diễn đ t trên, dễ dàng thấy rằng:  1.96  n≥  pˆ (1 − pˆ )  m  2 Do đó, số l ợng cỡ mẫu tùy thuộc vào độ sai số m và tỉ lệ p mà chúng ta muốn ớc tính. Độ sai số càng thấp, số l ợng c mẫu càng cao. 111 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Ví dụ 24: Chúng ta muốn ớc tính tỉ lệ đàn ông hút thuốc Việt Nam, sao cho ớc số không cao hơn hay thấp hơn 2% so với tỉ lệ thật trong toàn dân số. Một nghiên cứu tr ớc cho thấy tỉ lệ hút thuốc trong đàn ông ng i Việt có thể lên đến 70%. Câu hỏi đặt ra là chúng ta cần nghiên cứu trên bao nhiêu đàn ông để đ t yêu cầu trên. Trong ví dụ này, chúng ta có sai số m = 0.02, p̂ = 0.70, và số l ợng cỡ mẫu cần thiết cho nghiên cứu là:  1.96  n≥  0.7 × 0.3  0.02  2 Nói cách khác, chúng ta cần nghiên cứu ít nhất là 2017. Nếu chúng ta muốn gi m sai số từ 2% xuống 1% (tức m = 0.01) thì số l ợng đối t ợng sẽ là 8067! Chỉ cần thêm độ chính xác 1%, số l ợng mẫu có thể thêm hơn 6000 ng i. Do đó, vấn đề ớc tính cỡ mẫu ph i rất thận trọng, xem xét cân bằng giữa độ chính xác thông tin cần thu thập và chi phí. R không có hàm cho ớc tính cỡ mẫu cho một tỉ lệ, nh ng với công thức trên, b n đọc có thể viết một hàm để tính rất dễ dàng. 13.4.5 ớc tính cỡ mẫu cho so sánh hai tỉ lệ Nhiều nghiên cứu mang tính suy luận th ng có hai [hay nhiều hơn hai] nhóm để so sánh. Trong phần 15.4.2 chúng ta đã làm quen với ph ơng pháp ớc tính cỡ mẫu để so sánh hai số trung bình bằng kiểm định t. Đó là những ng i cứu mà tiêu chí là những biến số liên tục. Nh ng có nghiên cứu biến số không liên tục mà mang tính nhị phân nh tôi vừa bàn trong phần 15.4.3. Để so sánh hai tỉ lệ, ph ơng pháp kiểm định thông dụng nhất là kiểm định nhị phân (binomial test) hay Chi bình ph ơng (χ2 test). Trong phần này, tôi sẽ bàn qua cách tính cỡ mẫu cho hai lo i kiểm định thống kê này. Gọi hai tỉ lệ [mà chúng ta không biết nh ng muốn tìm hiểu] là p1 và p2 , và gọi ∆ = p1 – p2 . Gi thiết mà chúng ta muốn kiểm định là ∆ = 0. Lí thuyết đằng sau để ớc tính cỡ mẫu cho kiểm định gi thiết này khá r m rà, nh ng có thể tóm gọn bằng công thức sau đây: n= ( zα / 2 2 p (1 − p ) + zβ p1 (1 − p1 ) + p2 (1 − p2 ) ) 2 ∆ là trị số z của phân phối chuẩn cho xác suất α/2 (chẳng 2 Trong đó, p = ( p1 + p2 )/2, zα / 2 h n nh khi α = 0.05, thì zα / 2 = 1.96; khi α = 0.01, thì zα / 2 = 2.57), và zβ là trị số z của 112 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn phân phối chuẩn cho xác suất β (chẳng h n nh khi β = 0.10, thì zβ = 1.28; khi β = 0.20, thì zβ = 0.84). Ví dụ 25: Một thử nghiệm lâm sàng đối chứng ngẫu nhiên đ ợc thiết kế để đánh giá hiệu qu của một lo i thuốc chống gãy x ơng sống. Hai nhóm bệnh nhân sẽ đ ợc tuyển. Nhóm 1 đ ợc điều trị bằng thuốc, và nhóm 2 là nhóm đối chứng (không đ ợc điều trị). Các nhà nghiên cứu gi thiết rằng tỉ lệ gãy x ơng trong nhóm 2 là kho ng 10%, và thuốc có thể làm gi m tỉ lệ này xuống kho ng 6%. Nếu các nhà nghiên cứu muốn thử nghiệm gi thiết này với sai sót I là α = 0.01 và power = 0.90, bao nhiêu bệnh nhân cần ph i đ ợc tuyển mộ cho nghiên cứu? đây, chúng ta có ∆ = 0.10 – 0.06 = 0.04, và p = (0.10 + 0.06)/2 = 0.08. Với α = 0.01, zα / 2 = 2.57 và với power = 0.90, zβ = 1.28. Do đó, số l ợng bệnh nhân cần thiết cho mỗi nhóm là: ( 2.57 n= 2 × 0.08 × 0.92 + 1.28 0.1× 0.90 + 0.06 × 0.94 ( 0.04 ) 2 ) 2 = 1361 Nh vậy, công trình nghiên cứu này cần ph i tuyển ít nhất là 2722 bệnh nhân để kiểm định gi thiết trên. Hàm power.prop.test R có thể ứng dụng để tính cỡ mẫu cho tr ng hợp trên. Hàm power.prop.test cần những thông tin nh power, sig.level, p1, và p2. Trong ví dụ trên, chúng ta có thể viết: > power.prop.test(p1=0.10, p2=0.06, power=0.90, sig.level=0.01) Two-sample comparison of proportions power calculation n p1 p2 sig.level power alternative = = = = = = 1366.430 0.1 0.06 0.01 0.9 two.sided NOTE: n is number in *each* group Chú ý kết qu từ R có phần chính xác hơn (1366 đối t ợng cho mỗi nhóm) vì R dùng nhiều số lẽ cho tính toán hơn là tính “thủ công”. Tr ớc khi r i ch ơng này, tôi muốn nhân cơ hội này để nhấn m nh một lần nữa, ớc tính cỡ mẫu cho nghiên cứu là một b ớc cực kì quan trọng trong việc thiết kế một nghiên cứu cho có ý nghĩa khoa học, vì nó có thể quyết định thành b i của nghiên cứu. Tr ớc khi ớc tính cỡ mẫu nhà nghiên cứu cần ph i biết tr ớc (hay ít ra là có vài gi thiết cụ thể) về vấn đề mình quan tâm. ớc tính cỡ mẫu cần một số thông số nh đề cập đến 113 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn trong phần đầu của ch ơng, và nếu các thông số này không có thì không thể ớc tính đ ợc. Trong tr ng hợp một nghiên cứu hoàn toàn mới, tức ch a ai từng làm tr ớc đó, có thể các thông số về độ nh h ng và độ dao động đo l ng sẽ không có, và nhà nghiên cứu cần ph i tiến hành một số mô phỏng (simulation) hay một nghiên cứu sơ kh i để có những thông số cần thiết. Cách ớc tính cỡ mẫu bằng mô phỏng là một lĩnh vực nghiên cứu khá chuyên sâu, không nằm trong đề tài của sách này, nh ng b n đọc có thể tìm hiểu thêm ph ơng pháp này trong các sách giáo khoa về thống kê học cấp cao hơn. Trên đây là vài h ớng dẫn nhanh để b n đọc có thể sử dụng R cho phân tích số liệu và t o biểu đồ. Bài viết này thực chất là tóm l ợc từ cuốn Phân tích số liệu và tạo biểu đồ bằng R: hướng dẫn và thực hành, do Nhà xuất b n Đ i học Quốc gia Thành phố Hồ Chí Minh ấn hành vào năm 2006. Chi tiết về lí thuyết và một số ph ơng pháp khác nh phân tích sự kiện, xây dựng mô hình thống kê, mô phỏng, lập ch ơng, v.v… có thể tìm trong sách trên. 114 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 14. Tài liệu tham khảo Hiện nay, th viện sách về R còn t ơng đối khiêm tốn so với th viện cho các phần mềm th ơng m i nh SAS và SPSS. Tuy nhiên, trong th i đ i tiến bộ phi th ng về thông tin internet và toàn cầu hóa nh hiện nay, sách in và sách xuất b n trên website không còn là những khác nhau bao xa. Phần lớn chỉ dẫn về cách sử dụng R có thể tìm thấy r i rác đây đó trên các website từ các tr ng đ i học và website cá nhân trên khắp thế giới. Trong phần này tôi chỉ liệt kê một số sách mà b n đọc, nếu cần tham kh o thêm, nên tìm đọc. Trong quá trình viết cuốn sách mà b n đọc đang cầm trên tay, tôi cũng tham kh o một số sách và trang web mà tôi sẽ liệt kê sau đây với vài l i nhận xét cá nhân. Tài liệu tham kh o chính về R là bài báo của hai ng i sáng t o ra R: Ihaka R, Gentleman R. R: A language for data analysis and graphics. Journal of Computational and Graphical Statistics 1996; 5:299-314. • • • • • “Data Analysis and Graphics Using R – An Example Approach” (Nhà xuất b n Cambridge University Press, 2003) của John Maindonald nay đã xuất in l i lần thứ 2 với thêm một tác gi mới John Braun. Đây là cuốn sách rất có ích cho những ai muốn tìm hiểu và học về R. Năm ch ơng đầu của sách viết cho b n đọc ch a từng biết về R, còn các ch ơng sau thì viết cho các b n đọc đã biết cách sử dụng R thành th o. “Introductory Statistics With R” (Nhà xuất b n Springer, 2004) của Peter Dalgaard là một cuốn sách lo i căn b n cho R nhắm vào b n đọc ch a biết gì về R. Sách t ơng đối ngắn (chỉ kho ng 200 trang) nh ng khá đắt giá! “Linear Models with R” (Nhà xuất b n Chapman & Hall/CRC, 2004) của Julian Faraway. Sách hiện có thể t i từ internet xuống miễn phí t i website sau đây: http://www.stat.lsa.umich.edu/~faraway/book/pra.pdf hay http://cran.rproject.org/doc/contrib/Faraway-PRA.pdf. Tài liệu dài 213 trang. “R Graphics (Computer Science and Data Analysis)” (Nhà xuất b n Chapman & Hall/CRC, 2005) của Paul Murrell. Đây là cuốn sách chuyên về phân tích biểu đồ bằng R. Sách có rất nhiều mã để b n đọc có thể tự mình thiết kế các biểu đồ phức t p và … màu mè. “Modern Applied Statistics with S-Plus” (Nhà xuất b n Springer, 4th Edition, 2003) của W. N. Venables và B. D. Ripley đ ợc viết cho ngôn ngữ S-Plus nh ng tất c các lệnh và mã trong sách này đều có thể áp dụng cho R mà không cần thay đổi. (S-Plus là tiền thân của R, nh ng S-Plus là một phần mềm th ơng m i, còn R thì hoàn toàn miễn phí!) Đây là cuốn sách có thể nói là cuốn sách tham kh o cho tất c ai muốn phát triển thêm về R. Hai tác gi cũng là những chuyên gia có thẩm quyền về ngôn ngữ R. Sách dành cho b n đọc với trình độ cao về máy tính và thống kê học. 115 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn Các website quan trọng hay có ích về R • Rất nhiều tài liệu tham kh o có thể t i từ website chính thức của R sau đây: http://cran.R-project.org/other-docs.html Trong đó có một số tài liệu quan trọng nh “An Introduction to R” của W. N. Venables và B. D. Ripley. Địa chỉ internet: http://cran.r-project.org/doc/manuals/R-intro.pdf. • Vài tài liệu h ớng dẫn cách sử dụng R có thể t i (miễn phí) và tham kh o nh sau: “R for Beginners” (57 trang) của Emmanuel Paradis. Tài liệu đ ợc so n cho b n đọc mới làm quen với R. Địa chỉ internet: http://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf. “Using R for Data Analysis and Graphics: Introduction, Code and Commentary” (35 trang) của John Maindonald là một tóm l ợc các lệnh và hàm căn b n của R cho phân tích số liệu và biểu đồ. Chủ đề của tài liệu này rất gần với cuốn sách mà b n đang đọc. Địa chỉ internet: http://cran.r-project.org/doc/contrib/usingR.pdf “Statistical Analysis with R – a quick start” (46 trang) của Oleg Nenadic và Walter Zucchini. Web. Tài liệu h ớng dẫn cách ứng dụng R cho phân tích thống kê và biểu đồ. Địa chỉ internet: http://www.statoek.wiso.uni-goettingen.de/mitarbeiter/ogi/pub/r_workshop.pdf “A Brief Guide to R for Beginners in Econometrics” (31 trang) của M. Arai. Tài liệu chủ yếu so n cho giới phân tích thống kê kinh tế. Địa chỉ internet: http://people.su.se/~ma/R_intro “Notes on the use of R for psychology experiments and questionnaires” (39 trag) của Jonathan Baron và Yuelin Li. Web. Tài liệu đ ợc so n cho giới nghiên cứu tâm lí học và xã hội học. Có ví dụ về log-linear model và một số mô hình phân tích ph ơng sai trong tâm lí học. Địa chỉ internet: http://www.psych.upenn.edu/~baron/rpsych/rpsych.html • • StatsRus gồm một s u tập về các mẹo để sử dụng R hữu hiệu hơn (dài kho ng 80 trang). Địa chỉ internet: http://lark.cc.ukans.edu/pauljohn/R/statsRus.html Và sau cùng là một t i liệu “Hướng dẫn sử dụng R cho phân tích số liệu và biểu đồ” (kho ng 50 trang – th ng xuyên cập nhật hóa) do chính tôi viết bằng tiếng Việt. Website: www.R.ykhoa.net thực chất là tóm l ợc một số ch ơng chính của cuốn sách này. Trang web này còn có tất c các dữ liệu (datasets) và các mã sử trong trong sách để b n đọc có thể t i xuống máy tính cá nhân để sử dụng. 116 Phân tích số liệu và biểu đồ bằng R Nguyễn Văn Tuấn 15. Thuật ngữ dùng trong sách Tiếng Anh 95% confidence interval Akaike Information criterion (AIC) Analysis of covariance Analysis of variance (ANOVA) Bar chart Binomial distribution Box plot Categorical variable Clock chart Coefficient of correlation Coefficient of determination Coefficient of heterogeneity Combination Continuous variable Correlation Covariance Cross-over experiment Cumulative probability distribution Degree of freedom Determinant Discrete variable Dot chart Estimate Estimator Factorial analysis of variance Fixed effects Frequency Function Heterogeneity Histogram Homogeneity Hypothesis test Inverse matrix Latin square experiment Least squares method Linear Logistic regression analysis Linear regression analysis Tiếng Việt Kho ng tin cậy 95% Tiêu chuẩn thông tin Akaike Phân tích hiệp biến Phân tích ph ơng sai Biểu đồ thanh Phân phối nhị phân Biểu đồ hình hộp Biến thứ bậc Biểu đồ đồng hồ Hệ số t ơng quan Hệ số xác định bội Hệ số bất đồng nhất Tổ hợp Biến liên tục T ơng quan Hợp biến Thí nghiệm giao chéo Hàm phân phối tích lũy Bậc tự do Định thức Biến r i r c Biểu đồ điểm ớc số Hàm ớc l ợng thống kê Phân tích ph ơng sai cho thí nghiệm giai thừa nh h ng bất biến Tần số Hàm Bất đồng nhất Biểu đồ tần số Đồng nhất Kiểm định gi thiết Ma trận nghịch đ o Thí nghiệm hình vuông Latin Ph ơng pháp bình ph ơng nhỏ nhất Phân tích hồi qui tuyến tính logistic Phân tích hồi qui tuyến tính 117 Phân tích số liệu và biểu đồ bằng R Matrix Maximum likelihood method Mean Median Meta-analysis Missing value Model Multiple linear regression analysis Normal distribution Object Parameter Permutation Pie chart Poisson distribution Polynomial regression Probability Probability density distribution P-value Quantile Random effects Random variable Relative risk Repeated measure experiment Residual Residual mean square Residual sum of squares Scalar matrix Scatter plot Significance Simulation Standard deviation Standard error Standardized normal distribution Survival analysis Traposed matrix Variable Variance Weight Weighted mean Nguyễn Văn Tuấn Ma trận Ph ơng pháp hợp lí cực đ i Số trung bình Số trung vị Phân tích tổng hợp Giá trị không Mô hình Phân tích hồi qui tuyến tính đa biến Phân phối chuẩn Đối t ợng Thông số Hoán vị Biểu đồ hình tròn Phân phối Poisson Hồi qui đa thức Xác suất Hàm mật độ xác suất Trị số P Hàm định bậc nh h ng ngẫu nhiên Biến ngẫu nhiên Tỉ số nguy cơ t ơng đối Thí nghiệm tái đo l ng Phần d Trung bình bình ph ơng phần d Tổng bình ph ơng phần d Ma trận vô h ớng Biểu đồ tán x Có ý nghĩa thống kê Mô phỏng Độ lệch chuẩn Sai số chuẩn Phân phối chuẩn chuẩn hóa Phân tích biến cố Ma trận chuyển vị Biến (biến số) Ph ơng sai Trọng số Trung bình trọng số 118