›

Giáo trình Phân tích số liệu và tạo biểu đồ bằng R - Nguyễn Văn Tuấn

pdf 317 trang huongle 3230
Download
Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Phân tích số liệu và tạo biểu đồ bằng R - Nguyễn Văn Tuấn", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfgiao_trinh_phan_tich_so_lieu_va_tao_bieu_do_bang_r_nguyen_va.pdf

Nội dung text: Giáo trình Phân tích số liệu và tạo biểu đồ bằng R - Nguyễn Văn Tuấn

  1. NGUYỄN VĂN TUẤN PHÂN TÍCH SỐ LIỆU và TẠO BIỂU ĐỒ bằng Hướng dẫnthực hànhHướng dẫnthực hành
  2. Phân tích số liệu và tạo biểu đồ bằng hướng dẫn thực hành Mục lục 1 Lời nói đầu 2 Giới thiệu ngôn ngữ R 2.1 R là gì ? 2.2 Tải và cài đặt R vào máy tính 2.3 Package cho các phân tích đặc biệt 2.4 Khởi động và ngưng chạy R 2.5 “Văn phạm” ngôn ngữ R 2.6 Cách đặt tên trong R 2.7 Hỗ trợ trong R 2.8 Môi trường vận hành 3 Nhập dữ liệu 3.1 Nhập số liệu trực tiếp: c() 3.2 Nhập số liệu trực tiếp: edit(data.frame()) 3.3 Nhập số liệu từ một textfile: read.table() 3.4 Nhập số liệu từ Excel: read.csv 3.5 Nhập số liệu từ SPSS: read.spss 3.6 Tìm thông tin cơ bản về dữ liệu 4 Biên tập dữ liệu 4.1 Kiểm tra số liệu trống không: na.omit() 4.2 Tách rời dữ liệu: subset 4.3 Chiết số liệu từ một data .frame 4.4 Nhập hai data.frame thành một: merge 4.5 Mã hóa số liệu (data coding) 4.5.1 Mã hoá bằng hàm replace 4.5.2 Đổi một biến liên tục thành biến rời rạc 4.6 Chia một biến liên tục thành nhóm: cut 4.7 Tập hợp số liệu bằng cut2 (Hmisc) 1
  3. 5 Sử R cho các phép tính đơn giản và ma trận 5.1 Tính toán đơn giản 5.2 Số liệu về ngày tháng 5.3 Tạo dãy số bằng seq, rep và gl 5.4 Sử dụng R cho các phép tính ma trận 5.4.1 Chiết phần tử từ ma trận 5.4.2 Tính toán với ma trận 6 Tính toán xác suất và mô phỏng (simulation) 6.1 Tính toán đơn giản 6.1.1 Phép hoán vị (permutation) 6.1.2 Tổ hợp (combination) 6.2 Biến số ngẫu nhiên và hàm phân phối 6.3 Các hàm phân phối xác suất (probability distribution function) 6.3.1 Hàm phân phối nhị phân (Binomial distribution) 6.3.2 Hàm phân phối Poisson (Poisson distribution) 6.3.3 Hàm phân phối chuẩn (Normal distribution) 6.3.4 Hàm phân phối chuẩn chuẩn hóa (Standardized Normal distribution) 6.3.5 Hàm phân phối t, F và χ2 6.4. Mô phỏng (simulation) 6.4.1 Mô phỏng phân phối nhị phân 6.4.2 Mô phỏng phân phối Poisson 6.4.3 Mô phỏng phân phối χ2, t, F, gamma, beta, Weibull, Cauchy 6.5 Chọn mẫu ngẫu nhiên (random sampling) 7 Kiểm định giả thiết thống kê và ý nghĩa trị số P 7.1 Trị số P 7.2 Giả thiết khoa học và phản nghiệm 7.3 Ý nghĩa của trị số P qua mô phỏng 7.4 Vấn đề logic của trị số P 7.5 Vấn để kiểm định nhiều giả thiết (multiple tests of hypothesis) 8 Phân tích số liệu bằng biểu đồ 8.1 Môi trường và thiết kế biểu đồ 8.1.1 Nhiều biểu đồ cho một cửa sổ (windows) 8.1.2 Đặt tên cho trục tung và trục hoành 8.1.3 Cho giới hạn của trục tung và trục hoành 8.1.4 Thể loại và đường biểu diễn 8.1.5 Màu sắc, khung, và kí hiệu 8.1.6 Ghi chú (legend) 8.17 Viết chữ trong biểu đồ 2
  4. 8.2 Số liệu cho phân tích biểu đồ 8.3 Biểu đồ cho một biến số rời rạc (discrete variable): barplot 8.4. Biểu đồ cho hai biến số rời rạc (discrete variable): barplot 8.5 Biểu đồ hình tròn 8.6 Biểu đồ cho một biến số liên tục: stripchart và hist 8.6.1 Stripchart 8.6.2 Histogram 8.6.3 Biểu đồ hộp (boxplot) 8.6.4 Biểu đồ thanh (barchart) 8.6.5 Biểu đồ điểm (dotchart) 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) 8.8 Phân tích Biểu đồ cho nhiều biến: pairs 8.9 Một số biểu đồ “đa năng” 8.9.1 Biểu đồ tán xạ và hình hộp 8.9.2 Biểu đồ tán xạ với kích thước biến thứ ba 8.9.3 Biểu đồ thanh và xác suất tích lũy 8.9.4 Biểu đồ hình đồng hồ (clock plot) 8.9.5 Biểu đồ với sai số chuẩn (standard error) 8.9.6 Biểu đồ vòng (contour plot) 8.9.10 Biểu đồ với kí hiệu toán 9 Phân tích thống kê mô tả 9.0 Khái niệm về tổng thể (population) và mẫu (sample) 9.1 Thống kê mô tả: summary 9.2 Kiểm định xem một biến có phải phân phối chuẩn 9.3 Thống kê mô tả theo từng nhóm 9.4 Kiểm định t (t.test) 9.4.1 Kiểm định t một mẫu 9.4.2 Kiểm định t hai mẫu 9.5 So sánh phương sai (var.test) 9.6 Kiểm định Wilcoxon cho hai mẫu (wilcox.test) 9.7 Kiểm định t cho các biến số theo cặp (paired t-test, t.test) 9.8 Kiểm định Wilcoxon cho các biến số theo cặp (wilcox.test) 9.9 Tần số (frequency) 9.10 Kiểm định tỉ lệ (proportion test, prop.test, binom.test) 9.11 So sánh hai tỉ lệ (prop.test, binom.test) 9.12 So sánh nhiều tỉ lệ (prop.test, chisq.test) 9.12.1 Kiểm định Chi bình phương 9.12.2 Kiểm định Fisher 3
  5. 10 Phân tích hồi qui tuyến tính (regression analysis) 10.1 Hệ số tương quan 10.1.1 Hệ số tương quan Pearson 10.1.2 Hệ số tương quan Spearman 10.1.3 Hệ số tương quan Kendall 10.2 Mô hình của hồi qui tuyến tính đơn giản 10.2.1 Vài dòng lí thuyết 10.2.2 Phân tích hồi qui tuyến tính đơn giản bằng R 10.2.3 Giả định của phân tích hồi qui tuyến tính 10.2.4 Mô hình tiên đoán 10.3 Mô hình hồi qui tuyến tính đa biến (multiple linear regression) 10.4 Phân tích hồi qui đa thức (Polynomial regression analysis) 10.5 Xây dựng mô hình tuyến tính từ nhiều biến 10.6 Xây dựng mô hình tuyến tính bằng Bayesian Model Average (BMA) 11 Phân tích phương sai (analysis of variance) 11.1 Phân tích phương sai đơn giản (one-way analysis of variance - ANOVA) 11.1.1 Mô hình phân tích phương sai 11.1.2 Phân tích phương sai đơn giản với R 11.2 So sánh nhiều nhóm (multiple comparisons) và điều chỉnh trị số p 11.2.1 So sánh nhiều nhóm bằng phương pháp Tukey 11.2.2 Phân tích bằng biểu đồ 11.3 Phân tích bằng phương pháp phi tham số 11.4 Phân tích phương sai hai chiều (two-way analysis of variance - ANOVA) 11.4.1 Phân tích phương sai hai chiều với R 11.5 Phân tích hiệp biến (analysis of covariance - ANCOVA) 11.5.1 Mô hình phân tích hiệp biến 11.5.2 Phân tích bằng R 11.6 Phân tích phương sai cho thí nghiệm giai thừa (factorial experiment) 11.7 Phân tích phương sai cho thí nghiệm hình vuông Latin (Latin square experiment) 11.8 Phân tích phương sai cho thí nghiệm giao chéo (cross-over experiment) 11.9 Phân tích phương sai cho thí nghiệm tái đo lường (repeated measure experiment) 12 Phân tích hồi qui logistic (logistic regression analysis) 12.1 Mô hình hồi qui logistic 4
  6. 12.2 Phân tích hồi qui logistic bằng R 12.3 Ước tính xác suất bằng R 12.4 Phân tích hồi qui logistic từ số liệu giản lược bằng R 12.5 Phân tích hồi qui logistic đa biến và chọn mô hình 12.6 Chọn mô hình hồi qui logistic bằng Bayesian Model Average 12.7 Số liệu dùng cho phân tích 13 Phân tích biến cố (survival analysis) 13.1 Mô hình phân tích số liệu mang tính thời gian 13.2 Ước tính Kaplan-Meier bằng R 13.3 So sánh hai hàm xác suất tích lũy: kiểm định log-rank (log- rank test) 13.4 Kiểm định log-rank bằng R 13.5 Mô hình Cox (hay Cox’s proportional hazards model) 13.6 Xây dựng mô hình Cox bằng Bayesian Model Average (BMA) 14 Phân tích tổng hợp (meta-analysis) 14.1 Nhu cầu cho phân tích tổng hợp 14.2 Ảnh hưởng ngẫu nhiên và ảnh hưởng bất biến (Fixed- effects và Random-effects) 14.3 Qui trình của một phân tích tổng hợp 14.4 Phân tích tổng hợp ảnh hưởng bất biến cho một tiêu chí liên tục (Fixed-effects meta-analysis for a continuous outcome) 14.4.1 Phân tích tổng hợp bằng tính toán “thủ công” 14.4.2 Phân tích tổng hợp bằng R 14.5 Phân tích tổng hợp ảnh hưởng bất biến cho một tiêu chí nhị phân (Fixed-effects meta-analysis for a dichotomous outcome) 14.5.1 Mô hình phân tích 14.5.2 Phân tích bằng R 15 Ước tính cỡ mẫu (estimation of sample size) 15.1 Khái niệm về “power” 15.2 Thử nghiệm giả thiết thống kê và chẩn đoán bệnh 15.3 Số liệu để ước tính cỡ mẫu 15.4 Ước tính cỡ mẫu 15.4.1 Ước tính cỡ mẫu cho một chỉ số trung bình 15.4.2 Ước tính cỡ mẫu cho so sánh hai số trung bình 15.4.3 Ước tính cỡ mẫu cho phân tích phương sai 15.4.4 Ước tính cỡ mẫu cho ước tính một tỉ lệ 15.4.5 Ước tính cỡ mẫu cho so sánh hai tỉ lệ 16 Phụ lục 1: Lập trình và viết hàm bằng ngôn ngữ R 5
  7. 17 Phụ lục 2: Một số lệnh thông dụng trong R 18 Phụ lục 3: Thuật ngữ dùng trong sách 19 Lời bạt (tài liệu tham khảo và đọc thêm) 6
  8. CHƯƠNG I LỜI NÓI ĐẦU
  9. 1 Lời nói đầu Trái với quan điểm của nhiều người, thống kê là một bộ môn khoa học: Khoa học thống kê (Statistical Science). Các phương pháp phân tích dù dựa vào nền tảng của toán học và xác suất, nhưng đó chỉ là phần “kĩ thuật”, phần quan trọng hơn là thiết kế nghiên cứu và diễn dịch ý nghĩa dữ liệu. Người làm thống kê, do đó, không chỉ là người đơn thuần làm phân tích dữ liệu, mà phải là một nhà khoa học, một nhà suy nghĩ (“thinker”) về nghiên cứu khoa học. Chính vì thế, mà khoa học thống kê đóng một vai trò cực kì quan trọng, một vai trò không thể thiếu được trong các công trình nghiên cứu khoa học, nhất là khoa học thực nghiệm. Có thể nói rằng ngày nay, nếu không có thống kê thì các thử nghiệm gen với triệu triệu số liệu chỉ là những con số vô hồn, vô nghĩa. Một công trình nghiên cứu khoa học, cho dù có tốn kém và quan trọng cỡ nào, nếu không được phân tích đúng phương pháp sẽ không có ý nghĩa khoa học gì cả. Chính vì thế thế mà ngày nay, chỉ cần nhìn qua tất cả các tập san nghiên cứu khoa học trên thế giới, hầu như bất cứ bài báo y học nào cũng có phần “Statistical Analysis” (Phân tích thống kê), nơi mà tác giả phải mô tả cẩn thận phương pháp phân tích, tính toán như thế nào, và giải thích ngắn gọn tại sao sử dụng những phương pháp đó để hàm ý “bảo kê” hay tăng trọng lượng khoa học cho những phát biểu trong bài báo. Các tạp san y học có uy tín càng cao yêu cầu về phân tích thống kê càng nặng. Xin nhắc lại để nhấn mạnh: không có phần phân tích thống kê, bài báo không có ý nghĩa khoa học. Một trong những phát triển quan trọng nhất trong khoa học thống kê là ứng dụng máy tính cho phân tích và tính toán thống kê. Có thể nói không ngoa rằng không có máy tính, khoa học thống kê vẫn chỉ là một khoa học buồn tẻ khô khan, với những công thức rắc rối mà thiếu tính ứng dụng vào thực tế. Máy tính đã giúp khoa học thống kê làm một cuộc cách mạng lớn nhất trong lịch sử của bộ môn: đó là đưa khoa học thống kê vào thực tế, giải quyết các vấn đề gai góc nhất và góp phần làm phát triển khoa học thực nghiệm. Người viết còn nhớ hơn 20 năm về trước khi còn là một sinh viên theo học chương trình thạc sĩ thống kê ở Úc, một vị giáo sư khả kính kể một câu chuyện về nhà thống kê danh tiếng người Mĩ, Fred Mosteller, nhận được một hợp đồng nghiên cứu từ Bộ Quốc phòng Mĩ để cải tiến độ chính xác của vũ khí Mĩ vào thời Thế chiến thứ II, mà trong đó ông phải giải một bài toán thống kê gồm khoảng 30 thông số. Ông phải mướn 20 sinh viên sau đại học làm việc này: 10 sinh viên chỉ việc suốt ngày tính toán bằng tay; còn 10 sinh viên khác kiểm tra lại tính toán của 10 sinh viên kia. Công việc kéo dài gần một tháng trời. Ngày nay, với một máy tính cá nhân (personal computer) khiêm tốn, phân tích thống kê đó có thể giải trong vòng trên dưới 1 giây.
  10. Nhưng nếu máy tính mà không có phần mềm thì máy tính cũng chỉ là một đống sắt hay silicon “vô hồn” và vô dụng. Một phần mềm đã, đang và sẽ làm cách mạng thống kê là R. Phần mềm này được một số nhà nghiên cứu thống kê và khoa học trên thế giới phát triển và hoàn thiện trong khoảng 10 năm qua để sử dụng cho việc học tập, giảng dạy và nghiên cứu. Cuốn sách này sẽ giới thiệu bạn đọc cách sử dụng R cho phân tích thống kê và đồ thị. Tại sao R? Trước đây, các phần mềm dùng cho phân tích thống kê đã được phát triển và khá thông dụng. Những phần mềm nổi tiếng từ thời “xa xưa” như MINITAB, BMD-P đến những phần mềm tương đối mới như STATISTICA, SPSS, SAS, STAT, v.v thường rất đắt tiền (giá cho một đại học có khi lên đến hàng trăm ngàn đô-la hàng năm), một cá nhân hay thậm chí cho một đại học không khả năng mua. Nhưng R đã thay đổi tình trạng này, vì R hoàn toàn miễn phí. Trái với cảm nhận thông thường, miễn phí không có nghĩa là chất lượng kém. Thật vậy, chẳng những hoàn toàn miễn phí, R còn có khả năng làm tất cả (xin nói lại: tất cả), thậm chí còn hơn cả, những phân tích mà các phần mềm thương mại làm. R có thể tải xuống máy tính cá nhân của bất cứ cá nhân nào, bất cứ lúc nào, và bất cứ ở đâu trên thế giới. Chỉ vài phút cài đặt là R có thể đưa vào sử dụng. Chính vì thế mà đại đa số các đại học Tây phương và thế giới càng ngày càng chuyển sang sử dụng R cho học tập, nghiên cứu và giảng dạy. Trong xu hướng đó, cuốn sách này có một mục tiêu khiêm tốn là giới thiệu đến bạn đọc trong nước để kịp thời cập nhật hóa những phát triển về tính toán và phân tích thống kê trên thế giới. Cuốn sách này được soạn chủ yếu cho sinh viên đại học và các nhà nghiên cứu khoa học, những người cần một phần mềm để học thống kê, để phân tích số liệu, hay vẽ đồ thị từ số liệu khoa học. Cuốn sách này không phải là sách giáo khoa về lí thuyết thống kê, hay nhằm chỉ bạn đọc cách làm phân tích thống kê, nhưng sẽ giúp bạn đọc làm phân tích thống kê hữu hiệu hơn và hào hứng hơn. Mục đích chính của tôi là cung cấp cho bạn đọc những kiến thức cơ bản về thống kê, và cách ứng dụng R cho giải quyết vấn đề, và qua đó làm nền tảng để bạn đọc tìm hiểu hay phát triển thêm R. Tôi cho rằng, cũng như bất cứ ngành nghề nào, cách học phân tích thống kê hay nhất là tự mình làm phân tích. Vì thế, sách này được viết với rất nhiều ví dụ và dữ liệu thực. Bạn đọc có thể vừa đọc sách, vừa làm theo những chỉ dẫn trong sách (bằng cách gõ các lệnh vào máy tính) và sẽ thấy hào hứng hơn. Nếu bạn đọc đã có sẵn một dữ liệu nghiên cứu của chính mình thì việc học tập sẽ hữu hiệu hơn bằng cách ứng dụng ngay những phép tính trong sách. Đối với sinh viên, nếu chưa có số liệu sẵn, các bạn có thể dùng các phương pháp mô phỏng (simulation) để hiểu thống kê hơn. Khoa học thống kê ở nước ta tương đối còn mới, cho nên một số thuật ngữ chưa được diễn dịch một cách thống nhất và hoàn chỉnh. Vì thế, bạn đọc sẽ thấy đây đó trong sách một vài thuật ngữ “lạ”, và trong trường hợp này, tôi cố gắng kèm theo thuật ngữ gốc
  11. tiếng Anh để bạn đọc tham khảo. Ngoài ra, trong phần cuối của sách, tôi có liệt kê các thuật ngữ Anh – Việt đã được đề cập đến trong sách. Tất cả các dữ liệu sử dụng trong sách này đều có thể tải từ internet xuống máy tính cá nhân, hay có thể truy nhập trực tiếp qua trang web: Tôi hi vọng bạn đọc sẽ tìm thấy trong sách một vài thông tin bổ ích, một vài kĩ thuật hay phép tính có ích cho việc học tập, giảng dạy và nghiên cứu của mình. Nhưng có lẽ chẳng có cuốn sách nào hoàn thiện hay không có thiếu sót; thành ra, nếu bạn đọc phát hiện một sai sót trong sách, xin báo cho tôi biết qua điện thư t.nguyen@garvan.org.au hay rknguyen@gmail.com. Thành thật cám ơn các bạn đọc trước. Tôi muốn nhân dịp này cám ơn Tiến sĩ Nguyễn Hoàng Dzũng thuộc khoa Hóa, Đại học Bách khoa Thành phố Hồ Chí Minh, người đã gợi ý và giúp đỡ tôi in cuốn sách này ở trong nước. Tôi cám ơn Bác sĩ Nguyễn Đình Nguyên, ngừơi đã đọc một phần lớn bản thảo của cuốn sách, góp nhiều ý kiến thiết thực, và đã thiết kế bìa sách. Tôi cũng cám ơn Nhà xuất bản Đại học Bách khoa Thành phố Hồ Chí Minh đã giúp tôi in cuốn sách này. Bây giờ, tôi mời bạn đọc cùng đi với tôi một “hành trình thống kê” ngắn bằng R. Sydney, 31 Tháng Ba Năm 2006 Nguyễn Văn Tuấn
  12. CHƯƠNG II GIỚI THIỆU NGÔN NGỮ R
  13. 2 Giới thiệu ngôn ngữ R 2.1 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à đồ thị. 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. Hai người sáng tạo ra R là hai nhà thống kê học tên là Ross Ihaka và Robert Gentleman. Kể từ khi R ra đời, rất nhiều nhà nghiên cứu thống kê và toán học trên thế giới ủng hộ và tham gia vào việc phát triển R. Chủ trương của những người sáng tạo ra R là theo định hướng mở rộng (Open Access). Cũng một phần vì chủ trương này mà R hoàn toàn miễn phí. Bất cứ ai ở bất cứ nơi nào trên thế giới đều có thể truy nhập và tải toàn bộ mã nguồn của R về máy tính của mình để sử dụng. Cho đến nay, chỉ qua chưa đầy 5 năm phát triển, càng ngày càng có nhiều các 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 gần một triệu người sử dụng R, và con số này đang tăng theo cấp số nhân. Có thể nói trong vòng 10 năm nữa, chúng ta sẽ không cần đến các phần mềm thống kê đắt tiến như SAS, SPSS hay Stata (các phần mềm này rất đắt tiền, có thể lên đến 100.000 USD một năm) để phân tích thống kê nữa, vì tất cả các phân tích đó có thể tiến hành bằng R. 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. 2.2 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: 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à:
  14. R-2.2.1-win32.zip Tài liệu này khoảng 26 MB, và địa chỉ cụ thể để tải là: 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. 2.3 Package cho các phân tích đặc biệt 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ề máy tính và cài đặt. Địa chỉ để tải các package vẫn là: 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. Một số package cần tải về máy tính để sử dụng cho các ví dụ trong sách này là: Tên package Chức năng Trellis Dùng để vẽ đồ thị và làm cho đồ thị đẹp hơn lattice Dùng để vẽ đồ thị và làm cho đồ thị đẹp hơn Hmisc Một số phương pháp mô hình dữ liệu của F. Harrell Design Một số mô hình thiết kế nghiên cứu của F. Harrell Epi Dùng cho các phân tích dịch tễ học epitools Một package khác chuyên cho các phân tích dịch tễ học foreign Dùng để nhập dữ liệu từ các phần mềm khác như SPSS, Stata, SAS, v.v Rmeta Dùng cho phân tích tổng hợp (meta-analysis) meta Một package khác cho phân tích tổng hợp survival Chuyên dùng cho phân tích theo mô hình Cox (Cox’s proportional hazard model)
  15. splines Package cho survival vận hành Zelig Package dùng cho các phân tích thống kê trong lĩnh vực xã hội học genetics Package dùng cho phân tích số liệu di truyền học BMA Bayesian Model Average leaps Package dùng cho BMA 2.4 Khởi động và ngưng chạy R 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: R thường được sử dụng dưới dạng "command line", có nghĩa là chúng ta phải trực tiếp gõ lệnh vào cái prompt màu đỏ trên. Các lệnh phải tuân thủ nghiêm ngặt theo “văn phạm” và ngôn ngữ của R. Có thể nói toàn bộ bài viết này là nhằm hướng dẫn bạn đọc hiểu và viết theo ngôn ngữ của R. Một trong những văn phạm này là R phân biệt giữa Library và library. Nói cách khác, R phân biệt lệnh viết bằng chữ hoa hay chữ thường. Một văn phạm khác nữa là khi có hai chữ rời nhau, R thường dùng dấu chấm để
  16. thay vào khoảng trống, chẳng hạn như data.frame, t.test, read.table, v.v Điều này rất quan trọng, nếu không để ý sẽ làm mất thì giờ của người sử dụng. Nếu lệnh gõ ra đúng “văn phạm” thì R sẽ cho chúng ta một cái prompt khác hay cho ra kết quả nào đó (tùy theo lệnh); nếu lệnh không đúng văn phạm thì R sẽ cho ra một thông báo ngắn là không đúng hay không hiểu. Ví dụ, nếu chúng ta gõ: > x thì R sẽ hiểu và làm theo lệnh đó, rồi cho chúng ta một prompt khác: >. Nhưng nếu chúng ta gõ: > R is great R sẽ không “đồng ý” với lệnh này, vì ngôn ngữ này không có trong thư viện của R, một thông báo sau đây sẽ xuất hiện: Error: syntax error > Khi muốn rời khỏi R, chúng ta có thể đơn giản nhấn nút chéo (x) bên góc trái của window, hay gõ lệnh q(). 2.5 “Văn phạm” ngôn ngữ R “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. Chẳng hạn như: > reg 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, )
  17. 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 x y x lớn hơn y z = 1 p lớn hơn hoặc bằng 1 is.na(x) Có phải x là biến số missing A & B A và B (AND) A | B A hoặc B (OR) ! Không là (NOT) 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 myobject my object <- rnorm(10) Error: syntax error in "my object"
  18. 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 My.object.u my.object.L 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 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
  19. đâ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. Thay vì chọn mục trên, bạn đọc cũng có thể đơn giản lệnh: > help.start() và một cửa sổ sẽ xuất hiện chỉ dẫn toàn bộ hệ thống R. Hàm apropos cũng rất có ích vì nó cung cấp cho chúng ta tất cả các hàm trong R bắt đầu bằng kí tự mà chúng ta muốn tìm. Chẳng hạn như chúng ta muốn biết hàm nào trong R có kí tự “lm” thì chỉ đơn giản lệnh: > apropos(lm) Và R sẽ báo cáo các hàm với kí tự lm như sau có sẵn trong R: [1] ".__C__anova.glm" ".__C__anova.glm.null" ".__C__glm" [4] ".__C__glm.null" ".__C__lm" ".__C__mlm" [7] "anova.glm" "anova.glmlist" "anova.lm" [10] "anova.lmlist" "anova.mlm" "anovalist.lm" [13] "contr.helmert" "glm" "glm.control" [16] "glm.fit" "glm.fit.null" "hatvalues.lm" [19] "KalmanForecast" "KalmanLike" "KalmanRun" [22] "KalmanSmooth" "lm" "lm.fit" [25] "lm.fit.null" "lm.influence" "lm.wfit" [28] "lm.wfit.null" "model.frame.glm" "model.frame.lm"
  20. [31] "model.matrix.lm" "nlm" "nlminb" [34] "plot.lm" "plot.mlm" "predict.glm" [37] "predict.lm" "predict.mlm" "print.glm" [40] "print.lm" "residuals.glm" "residuals.lm" [43] "rstandard.glm" "rstandard.lm" "rstudent.glm" [46] "rstudent.lm" "summary.glm" "summary.lm" [49] "summary.mlm" "kappa.lm" 2.8 Môi trường vận hành Dữ liệu phải được chứa trong một khu vực (directory) của máy tính. Trước khi sử dụng R, có lẽ cách hay nhất là tạo ra một directory để chứa dữ liệu, chẳng hạn như c:\works\stats. Để R biết dữ liệu nằm ở đâu, chúng ta sử dụng lệnh setwd (set working directory) như sau: > setwd(“c:/works/stats”) Lệnh trên báo cho R biết là dữ liệu sẽ chứa trong directory có tên là c:\works\stats. Chú ý rằng, R dùng forward slash “/” chứ không phải backward slash “\” như trong hệ thống Windows. Để biết hiện nay, R đang “làm việc” ở directory nào, chúng ta chỉ cần lệnh: > getwd() [1] "C:/Program Files/R/R-2.2.1" Cái prompt mặc định của R là “>”. Nhưng nếu chúng ta muốn có một prompt khác theo cá tính cá nhân, chúng ta có thể thay thế dễ dàng: > options(prompt=”R> ”) R> Hay: > options(prompt="Tuan> ") Tuan> Màn ảnh R mặc định là 80 characters, nhưng nếu chúng ta muốn màn ảnh rộng hơn, thì chỉ cần ra lệnh: > options(width=100) Hay muốn R trình bày các số liệu ở dạng 3 số thập phân: > options(scipen=3)
  21. Các lựa chọn và thay đổi này có thể dùng lệnh options(). Để biết các thông số hiện tại của R là gì, chúng ta chỉ cần lệnh: > options() Tìm hiểu ngày tháng: > Sys.Date() [1] "2006-03-31" Nếu bạn đọc cần thêm thông tin, một số tài liệu trên mạng (viết bằng tiếng Anh) cũng rất có ích. Các tài liệu này có thể tải xuống máy miễn phí: R for beginners (của Emmanuel Paradis): Using R for data analysis and graphics (của John Maindonald):
  22. CHƯƠNG III NHẬP DỮ LIỆU
  23. 3 Nhập dữ liệu 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: 3.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 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 Chúng ta có thể sử dụng function có tên c như sau: > age insulin <- c(16.5,10.8,32.3,19.3,14.2,11.3,15.5,15.8,16.2,11.2) 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ờ
  24. 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 tuan Và R sẽ báo cáo: age insulin 1 50 16.5 2 62 10.8 3 60 32.3 4 40 19.3 5 48 14.2 6 47 11.3 7 57 15.5 8 70 15.8 9 48 16.2 10 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\stats”, chúng ta cần gõ như sau: > setwd(“c:/works/stats”) > save(tuan, file=”tuan.rda”) 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\stats”. 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 đó. 3.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,
  25. 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 window 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 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. 3.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\stats. 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 sex age bmi hdl ldl tc tg 1 Nam 57 17 5.000 2.0 4.0 1.1 2 Nu 64 18 4.380 3.0 3.5 2.1
  26. 3 Nu 60 18 3.360 3.0 4.7 0.8 4 Nam 65 18 5.920 4.0 7.7 1.1 5 Nam 47 18 6.250 2.1 5.0 2.1 6 Nu 65 18 4.150 3.0 4.2 1.5 7 Nam 76 19 0.737 3.0 5.9 2.6 8 Nam 61 19 7.170 3.0 6.1 1.5 9 Nam 59 19 6.942 3.0 5.9 5.4 10 Nu 57 19 5.000 2.0 4.0 1.9 11 Nu 63 20 4.217 5.0 6.2 1.7 12 Nam 51 20 4.823 1.3 4.1 1.0 13 Nu 60 20 3.750 1.2 3.0 1.6 14 Nam 42 20 1.904 0.7 4.0 1.1 15 Nam 64 20 6.900 4.0 6.9 1.5 16 Nu 49 20 0.633 4.1 5.7 1.0 17 Nu 44 21 5.530 4.3 5.7 2.7 18 Nu 45 21 6.625 4.0 5.3 3.9 19 Nu 80 21 5.960 4.3 7.1 3.0 20 Nu 48 21 3.800 4.0 3.8 3.1 21 Nu 61 21 5.375 3.1 4.3 2.2 22 Nu 45 21 3.360 3.0 4.8 2.7 23 Nu 70 21 5.000 1.7 4.0 1.1 24 Nu 51 21 2.608 2.0 3.0 0.7 25 Nam 63 22 4.130 2.1 3.1 1.0 26 Nam 54 22 5.000 4.0 5.3 1.7 27 Nu 57 22 6.235 4.1 5.3 2.9 28 Nam 70 22 3.600 4.0 5.4 2.5 29 Nu 47 22 5.625 4.2 4.5 6.2 30 Nu 60 22 5.360 4.2 5.9 1.3 31 Nu 60 22 6.580 4.4 5.6 3.3 32 Nam 50 22 7.545 4.3 8.3 3.0 33 Nam 60 22 6.440 2.3 5.8 1.0 34 Nu 55 22 6.170 6.0 7.6 1.4 35 Nu 74 23 5.270 3.0 5.8 2.5 36 Nam 48 23 3.220 3.0 3.1 0.7 37 Nu 46 23 5.400 2.6 5.4 2.4 38 Nam 49 23 6.300 4.4 6.3 2.4 39 Nu 69 23 9.110 4.3 8.2 1.4 40 Nu 72 23 7.750 4.0 6.2 2.7 41 Nam 51 23 6.200 3.0 6.2 2.4 42 Nu 58 23 7.050 4.1 6.7 3.3 43 nam 60 24 6.300 4.4 6.3 2.0 44 Nam 45 24 5.450 2.8 6.0 2.6 45 Nam 63 24 5.000 3.0 4.0 1.8 46 Nu 52 24 3.360 2.0 3.7 1.2 47 Nam 64 24 7.170 1.0 6.1 1.9 48 Nam 45 24 7.880 4.0 6.7 3.3 49 Nu 64 25 7.360 4.6 8.1 4.0 50 Nu 62 25 7.750 4.0 6.2 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/stats”)
  27. > chol chol Hay > names(chol) R sẽ cho biết có các cột như sau trong dữ liệu (name 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") 3.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 Age Sex Ethnicity IGFI IGFBP3 ALS PINP ICTP P3NP 1 18 1 1 148.27 5.14 316.00 61.84 5.81 4.21 2 28 1 1 114.50 5.23 296.42 98.64 4.96 5.33 3 20 1 1 109.82 4.33 269.82 93.26 7.74 4.56 4 21 1 1 112.13 4.38 247.96 101.59 6.66 4.61 5 28 1 1 102.86 4.04 240.04 58.77 4.62 4.95 6 23 1 4 129.59 4.16 266.95 48.93 5.32 3.82 7 20 1 1 142.50 3.85 300.86 135.62 8.78 6.75 8 20 1 1 118.69 3.44 277.46 79.51 7.19 5.11 9 20 1 1 197.69 4.12 335.23 57.25 6.21 4.44 10 20 1 1 163.69 3.96 306.83 74.03 4.95 4.84
  28. 11 22 1 1 144.81 3.63 295.46 68.26 4.54 3.70 12 27 0 2 141.60 3.48 231.20 56.78 4.47 4.07 13 26 1 1 161.80 4.10 244.80 75.75 6.27 5.26 14 33 1 1 89.20 2.82 177.20 48.57 3.58 3.68 15 34 1 3 161.80 3.80 243.60 50.68 3.52 3.35 16 32 1 1 148.50 3.72 234.80 83.98 4.85 3.80 17 28 1 1 157.70 3.98 224.80 60.42 4.89 4.09 18 18 0 2 222.90 3.98 281.40 74.17 6.43 5.84 19 26 0 2 186.70 4.64 340.80 38.05 5.12 5.77 20 27 1 2 167.56 3.56 321.12 30.18 4.78 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)” Sau khi xong, chúng ta sẽ có một file với tên “excel.csv” trong directory “c:\works\stats”. Việc thứ hai là vào R và ra những lệnh sau đây: > setwd(“c:/works/stats”) > 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.
  29. 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") 3.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\stats, 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: > library(foreign) Việc thứ hai là lệnh read.spss: > setwd(“c:/works/stats”) > testo save(testo, file="testo.rda") 3.6 Thông tin cơ bản 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
  30. 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" "tg" • 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 Nu 1 21 28 Kết quả cho thấy dữ liệu này có 21 nam và 28 nữ.
  31. CHƯƠNG IV BIÊN TẬP DỮ LIỆU
  32. 4 Biên tập dữ liệu Biên tập số liệu ở đây không có nghĩa là thay đổi số liệu gốc (vì đó là một tội lớn, một sự gian dối trong khoa học không thể chấp nhận được), mà chỉ có nghĩa tổ chức số liệu sao cho R có thể phân tích một cách hữu hiệu. Nhiều khi trong phân tích thống kê, chúng ta cần phải tập trung số liệu thành một nhóm, hay tách rời thành từng nhóm, hay thay thế từ kí tự (characters) sang số (numeric) cho tiện việc tính toán. Trong chương này, tôi sẽ bàn qua một số lệnh căn bản cho việc biên tập số liệu. 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/stats”) > chol attach(chol) 4.1 Kiểm tra số liệu trống không (missing value) Trong nghiên cứu, vì nhiều lí do số liệu không thể thu thập được cho tất cả đối tượng, hay không thể đo lường tất cả biến số cho một đối tượng. Trong trường hợp đó, số liệu trống được xem là “missing value” (mà tôi tạm dịch là số liệu trống không). R xem các số liệu trống không là NA. Có một số kiểm định thống kê đòi hỏi các số liệu trống không phải được loại ra (vì không thể tính toán được) trước khi phân tích. R có một lệnh rất có ích cho việc này: na.omit, và cách sử dụng như sau: > chol.new nam nu <- subset(chol, sex==”Nu”)
  33. 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 =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 =60 & sex==”Nam”) > dim(n60) [1] 9 8 4.3 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 data3 print(data3) id sex tc 1 1 Nam 4.0 2 2 Nu 3.5 3 3 Nu 4.7 4 4 Nam 7.7 5 5 Nam 5.0 6 6 Nu 4.2 7 7 Nam 5.9 8 8 Nam 6.1
  34. 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). 4.4 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 sex tg 1 Nam 1.1 2 Nu 2.1 3 Nu 0.8 4 Nam 1.1 5 Nam 2.1 6 Nu 1.5 7 Nam 2.6 8 Nam 1.5 9 Nam 5.4 10 Nu 1.9 11 Nu 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 d id sex.x tc sex.y tg
  35. 1 1 Nam 4.0 Nam 1.1 2 2 Nu 3.5 Nu 2.1 3 3 Nu 4.7 Nu 0.8 4 4 Nam 7.7 Nam 1.1 5 5 Nam 5.0 Nam 2.1 6 6 Nu 4.2 Nu 1.5 7 7 Nam 5.9 Nam 2.6 8 8 Nam 6.1 Nam 1.5 9 9 Nam 5.9 Nam 5.4 10 10 Nu 4.0 Nu 1.9 11 11 NA Nu 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”). 4.5 Mã hóa 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 diagnosis diagnosis[bmd diagnosis[bmd > -2.5 & bmd diagnosis[bmd > -1.0] data <- data.frame(bmd, diagnosis) # liệt kê để kiểm tra xem lệnh có hiệu quả không
  36. > data bmd diagnosis 1 -0.92 3 2 0.21 3 3 0.17 3 4 -3.21 1 5 -1.80 2 6 -2.60 1 7 -2.00 2 8 1.71 3 9 2.12 3 10 -2.11 2 4.5.1 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 -2.5 & bmd diagnosis -1.0, 3) 4.5.2 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 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:
  37. > 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ả. 4.6 Chia nhóm bằng cut Với một biến liên tục, chúng ta có thể chia thành nhiều nhóm bằng hàm cut. Ví dụ, chúng ta có biến age như sau: > age cut(age, 2) [1] (7.96,29.5] (7.96,29.5] (7.96,29.5] (29.5,51] (7.96,29.5] (7.96,29.5] (7.96,29.5] (7.96,29.5] [9] (7.96,29.5] (29.5,51] (7.96,29.5] (7.96,29.5] (7.96,29.5] (29.5,51] (29.5,51] Levels: (7.96,29.5] (29.5,51] cut chia biến age thành 2 nhóm: nhóm 1 tuổi từ 7.96 đến 29.5; nhóm 2 từ 29.5 đến 51. Chúng ta có thể đếm số đối tượng trong từng nhóm tuổi bằng hàm table như sau: > table(cut(age, 2)) (7.96,29.5] (29.5,51] 11 4 > ageg ageg table(ageg) ageg low medium high 10 2 3 Tất nhiên, chúng ta cũng có thể chia age thành 4 nhóm (quartiles) bằng cách cho những thông số 0, 0.25, 0.50 và 0.75 như sau: cut(age, breaks=quantiles(age, c(0, 0.25, 0.50, 0.75, 1)), labels=c(“q1”, “q2”, “q3”, “q4”),
  38. include.lowest=TRUE) cut(age, breaks=quantiles(c(0, 0.25, 0.50, 0.75, 1)), labels=c(“q1”, “q2”, “q3”, “q4”), include.lowest=TRUE) 4.7. Tập hợp số liệu bằng cut2 (Hmisc) Hàm cut trên chia biến số theo giá trị của biến, chứ không dựa vào số mẫu, cho nên số lượng mẫu trong từng nhóm không bằng nhau. Tuy nhiên, 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ố nhưng số mẫu bằng hay tương đương nhau. 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 với số mẫu 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 # chia biến số bmd thành 2 nhóm và để trong đối tượng group > group 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 table(group) group [-3.21,-1.80) [-1.80, 0.21) [ 0.21, 2.12] 4 3 3
  39. CHƯƠNG V TÍNH TOÁN ĐƠN GIẢN VÀ MA TRẬN
  40. 5 Dùng R cho các phép tính đơn giản và ma trậ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. 5.1 Tính toán đơn giản Cộng hai số hay nhiều số với nhau: Cộng và trừ: > 15+2997 > 15+2997-9768 [1] 3012 [1] -6756 Nhân và chia Số lũy thừa: (25 – 5)3 > -27*12/21 > (25 - 5)^3 [1] -15.42857 [1] 8000 Căn số bậc hai: 10 Số pi (π) > sqrt(10) > pi [1] 3.162278 [1] 3.141593 > 2+3*pi [1] 11.42478 Logarit: loge Logarit: log10 > log(10) > log10(100) [1] 2.302585 [1] 2 Số mũ: e2.7689 Hàm số lượng giác > exp(2.7689) > cos(pi) [1] 15.94109 [1] -1 > log10(2+3*pi) [1] 1.057848 Vector > exp(x/10) > x x 1.491825 1.822119 2.013753 1.822119 [1] 2 3 1 5 4 6 7 6 8 [9] 2.225541 > sum(x) > exp(cos(x/10)) [1] 2.664634 2.599545 2.704736 2.405 [1] 42 2.511954 2.282647 2.148655 2.282647 [9] 2.007132 > x*2
  41. [1] 4 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 2 2 2 2 n + 2 + 3 + 4 + 5 = ? 2 > x sum(x^2) i=1 [1] 55 > x 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 n 2 chuẩn (standard deviation): ()x − xn/ = ? n ∑ i 2 i=1 2 Phương sai: sxxn= ∑()()i −−/1= ? > x sum((x-mean(x))^2)/length(x) > x var(x) [1] 2.5 Trong công thức trên, length(x) có Độ lệch chuẩn: s2 : nghĩa là tổng số phần tử (elements) trong > sd(x) vector x. [1] 1.581139 5.2 Số liệu về ngày tháng Trong phân tích thống kê, các số liệu ngày tháng có khi là một vấn đề nan giải, vì có rất nhiều cách để mô tả các dữ liệu này. Chẳng hạn như 01/02/2003, có khi người ta viết 1/2/2003, 01/02/03, 01FEB2003, 2003-02-01, v.v Thật ra, có một qui luật chuẩn để viết số liệu ngày tháng là tiêu chuẩn ISO 8601 (nhưng rất ít ai tuân theo!) Theo qui luật này, chúng ta viết: 2003-02-01 Lí do đằng sau cách viết này là chúng ta viết số với đơn vị lớn nhất trước, rồi dần dần đến đơn vị nhỏ nhất. Chẳng hạn như với số “123” thì chúng ta biết ngay rằng “một trăm hai mươi ba”: bắt đầu là hàng trăm, rồi đến hàng chục, v.v Và đó cũng là cách viết ngày tháng chuẩn của R. > date1 date2 <- as.Date(“06/03/01”, format=”%y/%m/%d”) Chú ý chúng ta nhập hai số liệu khác nhau về thứ tự ngày tháng năm, nhưng chúng ta cũng cho biết cụ thể cách đọc bằng %d (ngày), %m (tháng), và %y (năm). Chúng ta có thể tính số ngày giữa hai thởi điểm:
  42. > days days Time difference of 28 days Chúng ta cũng có thể tạo một dãy số liệu ngày tháng như sau: > seq(as.Date(“2005-01-01”), as.Date(“2005-12-31”), by=”month”) [1] "2005-01-01" "2005-02-01" "2005-03-01" "2005-04-01" "2005-05-01" [6] "2005-06-01" "2005-07-01" "2005-08-01" "2005-09-01" "2005-10-01" [11] "2005-11-01" "2005-12-01" > seq(as.Date(“2005-01-01”), as.Date(“2005-12-31”), by=”2 weeks”) [1] "2005-01-01" "2005-01-15" "2005-01-29" "2005-02-12" "2005-02-26" [6] "2005-03-12" "2005-03-26" "2005-04-09" "2005-04-23" "2005-05-07" [11] "2005-05-21" "2005-06-04" "2005-06-18" "2005-07-02" "2005-07-16" [16] "2005-07-30" "2005-08-13" "2005-08-27" "2005-09-10" "2005-09-24" [21] "2005-10-08" "2005-10-22" "2005-11-05" "2005-11-19" "2005-12-03" [26] "2005-12-17" "2005-12-31" 5.3 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 x [1] 1 2 3 4 5 6 7 8 9 10 11 12 > seq(12) [1] 1 2 3 4 5 6 7 8 9 10 11 12 • Tạo ra một vector số từ 12 đến 5: > x x [1] 12 11 10 9 8 7 6 5 > seq(12,7) [1] 12 11 10 9 8 7 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:
  43. • Tạo ra một vector số từ 4 đến 6 với khoảng cách bằng 0.25: > seq(4, 6, 0.25) [1] 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 • 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(length=10, from=2, to=15) [1] 2.000000 3.444444 4.888889 6.333333 7.777778 9.222222 10.666667 12.111111 13.555556 15.000000 Á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: > rep(10, 3) [1] 10 10 10 • Tạo ra số 1 đến 4, 3 lần: > rep(c(1:4), 3) [1] 1 2 3 4 1 2 3 4 1 2 3 4 • Tạo ra số 1.2, 2.7, 4.8, 5 lần: > 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 • Tạo ra số 1.2, 2.7, 4.8, 5 lần: > 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)
  44. [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: > 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 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.4 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: 147  A = 258  369
  45. Và với R: > y A 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 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ụ: > y A 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 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 # cho các phần tử đường chéo bằng 1
  46. > diag(A) 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 5.4.1 Chiết phần tử từ ma trận > y A 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] [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
  47. > # 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 5.4.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 A [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12 > B B [,1] [,2] [,3] [,4] [1,] -1 -4 -7 -10 [2,] -2 -5 -8 -11 [3,] -3 -6 -9 -12 Chúng ta có thể cộng A+B: > C C [,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 0 0 0 0 [3,] 0 0 0 0 Hay A-B: > D 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:
  48. 147 123   A = 258 và B = 456   369 789 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 A B AB 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 BA [,1] [,2] [,3] [1,] 14 32 50 [2,] 32 77 122 [3,] 50 122 194 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: 34xx+= 4 12 xx12+=62 Hệ phương trình này có thể viết bằng kí hiệu ma trận: AX = Y, trong đó: 34 x1 4 A = , X = , và Y =  16 x2 2 Nghiệm của hệ phương trình này là: X = A-1Y, hay trong R: > A Y X X [,1] [1,] 1.1428571 [2,] 0.1428571
  49. 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 E [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > det(E) [1] 0 Nhưng ma trận F sau đây thì có thể đảo nghịch: > F 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
  50. 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ể. Trong vài chương sau, tôi sẽ quay lại vấn đề này chi tiết hơn. 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à: cùng với tài liệu chỉ dẫn cách sử dụng (dài khoảng 80 trang):
  51. CHƯƠNG VI TÍNH TOÁN XÁC SUẤT
  52. 6 Tính toán xác suất và mô phỏng (simulation) Xác suất là nền tảng của phân tích thống kê. Tất cả các phương pháp phân tích số liệu và suy luận thống kê đều dựa vào lí thuyết xác suất. Lí thuyết xác suất quan tâm đến việc mô tả và thể hiện qui luật phân phối của một biến số ngẫu nhiên. “Mô tả” ở đây trong thực tế cũng có nghĩa đơn giản là đếm những trường hợp hay khả năng xảy ra của một hay nhiều biến. Chẳng hạn như khi chúng ta chọn ngẫu nhiên 2 đối tượng , và nếu 2 đối tượng này có thể được phân loại bằng hai đặc tính như giới tính và sở thích, thì vấn đề đặt ra là có bao nhiêu tất cả “phối hợp” giữa hai đặc tính này. Hay đối với một biến số liên tục như huyết áp, mô tả có nghĩa là tính toán các chỉ số thống kê của biến như trị số trung bình, trung vị, phương sai, độ lệch chuẩn, v.v Từ những chỉ số mô tả, lí thuyết xác suất cung cấp cho chúng ta những mô hình để thiết lập các hàm phân phối cho các biến số đó. Trong chương này, tôi sẽ bàn qua hai lĩnh vực chính là phép đếm và các hàm phân phối. 6.1 Các phép đếm 6.1.1 Phép hoán vị (permutation). Theo định nghĩa, hoán vị n phần tử là cách sắp xếp n phần tử theo một thứ tự định sẵn. Định nghĩa này thật là khó hiểu, chẳng khác gì đố! Có lẽ một ví dụ cụ thể sẽ làm rõ định nghĩa hơn. Hãy tưởng tượng một trung tâm cấp cứu có 3 bác sĩ (x, y và z), và có 3 bệnh nhân (a, b và c) đang ngồi chờ được khám bệnh. Cả ba bác sĩ đều có thể khám bất cứ bệnh nhân a, b hay c. Câu hỏi đặt ra là có bao nhiêu cách sắp xếp bác sĩ – bệnh nhân? Để trả lời câu hỏi này, chúng ta xem xét vài trường hợp sau đây: • Bác sĩ x có 3 lựa chọn: khám bệnh nhân a, b hoặc c; • Khi bác sĩ x đã chọn một bệnh nhân rồi, thì bác sĩ y có hai lựa chọn còn lại; • Và sau cùng, khi 2 bác sĩ kia đã chọn, bác sĩ z chỉ còn 1 lựa chọn. • Tổng cộng, chúng ta có 6 lựa chọn. Một ví dụ khác, trong một buổi tiệc gồm 6 bạn, hỏi có bao nhiêu cách sắp xếp cách ngồi trong một bàn với 6 ghế? Qua cách lí giải của ví dụ trên, đáp số là: 6.5.4.3.2.1 = 720 cách. (Chú ý dấu “.” có nghĩa là dấu nhân hay tích số). Và đây chính là phép đếm hoán vị. 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à: nnnn!=−−()() 1 2( n −×× 3) 1. Trong R cách tính này rất đơn giản với lệnh prod() như sau: 1
  53. • Tìm 3! > prod(3:1) [1] 6 • Tìm 10! > prod(10:1) [1] 3628800 • 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 6.1.2 Tổ hợp (combination). Tổ hợp n phần tử chập k là mọi tập hợp con gồm k phần tử của tập hợp n phần tử. Định nghĩa này phải nói là rất khó hiểu và rườm rà! Cách dễ hiểu nhất là qua một ví dụ như sau: Cho 3 người (hãy cho là A, B, và C) ứng viên vào 2 chức chủ tịch và phó chủ tịch, hỏi: có bao nhiêu cách để chọn 2 chức này trong số 3 người đó. Chúng ta có thể tưởng tượng có 2 ghế mà phải chọn 3 người: Cách chọn Chủ tịch Phó chủ tịch 1 A B 2 B A 3 A C 4 C A 5 B C 6 C B Như vậy có 6 cách chọn. Nhưng chú ý rằng cách chọn 1 và 2 trong thực tế chỉ là 1 cặp, và chúng ta chỉ có thể đếm là 1 (chứ không 2 được). Tương tự, 3 và 4, 5 và 6 cũng chỉ có thể đếm là 1 cặp. Tổng cộng, chúng ta có 3 cách chọn 3 người cho 2 chức vụ. Đáp số này được gọi là tổ hợp. Thật ra tổng số lần chọn có thể tính bằng công thức sau đây: 3 3! 6 = ==3 lần. 2 2!() 3− 2 ! 2 Nói chung, số lần chọn k người từ n người là: n n! = k knk!!()− 2
  54. n n Công thức này cũng có khi viết là Ck thay vì . Với R, phép tính này rất đơn giản k bằng hàm choose(n, k). Sau đâ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ười được đắc cử vào hai chức vụ: > 1/choose(5, 2) [1] 0.1 6.2 Biến số ngẫu nhiên và hàm phân phối Phần lớn phân tích thống kê dựa vào các luật phân phối xác suất để suy luận. Hai chữ “phân phối” (distribution) có lẽ cũng cần vài dòng giải thích ở đây. Nếu chúng ta chọn ngẫu nhiên 10 bạn trong một lớp học và ghi nhận chiều cao và giới tính của 10 bạn đó, chúng ta có thể có một dãy số liệu như sau: 1 2 3 4 5 6 7 8 9 10 Giới tính Nữ Nữ Nam Nữ Nữ Nữ Nam Nam Nữ Nam Chiều cao (cm) 156 160 175 145 165 158 170 167 178 155 Nếu tính gộp chung lại, chúng ta có 6 bạn gái và 4 bạn trai. Nói theo phần trăm, chúng ta có 60% nữ và 40% nam. Nói theo ngôn ngữ xác suất, xác suất nữ là 0.6 và nam là 0.4. Về chiều cao, chúng ta có giá trị trung bình là 162.9 cm, với chiều cao thấp nhất là 155 cm và cao nhất là 178 cm. Nói theo ngôn ngữ thống kê xác suất, biến số giới tính và chiều cao là hai biến số ngẫu nhiên (random variable). Ngẫu nhiên là vì chúng ta không đoán trước một cách chính xác các giá trị này, nhưng chỉ có thể đoán giá trị tập trung, giá trị trung bình, và độ dao động của chúng. Biến giới tính chỉ có hai “giá trị” (nam hay nữ), và được gọi là biến không liên tục, hay biến rời rạc (discrete variable), hay biến thứ bậc (categorical variable). Còn biến chiều cao có thể có bất cứ giá trị nào từ thấp đến cao, và do đó có tên là biến liên tục (continuous variable). 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 3
  55. 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 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 Mật độ Tích lũy Định bậc Mô phỏng phối Chuẩn dnorm(x, mean, pnorm(q, mean, sd) qnorm(p, mean, sd) rnorm(n, mean, sd) sd) Nhị phân dbinom(k, n, p) pbinom(q, n, p) qbinom (p, n, p) rbinom(k, n, prob) Poisson dpois(k, lambda) ppois(q, lambda) qpois(p, lambda) rpois(n, lambda) Uniform dunif(x, min, punif(q, min, max) qunif(p, min, max) runif(n, min, max) max) Negative dnbinom(x, k, p) pnbinom(q, k, p) qnbinom (p,k,prob) rbinom(n, n, prob) binomial Beta dbeta(x, shape1, pbeta(q, shape1, qbeta(p, shape1, rbeta(n, shape1, shape2) shape2) shape2) shape2) Gamma dgamma(x, shape, gamma(q, shape, qgamma(p, shape, rgamma(n, shape, rate, scale) rate, scale) rate, scale) rate, scale) Geometric dgeom(x, p) pgeom(q, p) qgeom(p, prob) rgeom(n, prob) Exponential dexp(x, rate) pexp(q, rate) qexp(p, rate) rexp(n, rate) Weibull dnorm(x, mean, pnorm(q, mean, sd) qnorm(p, mean, sd) rnorm(n, mean, sd) sd) Cauchy dcauchy(x, pcauchy(q, qcauchy(p, rcauchy(n, location, scale) location, scale) location, scale) location, scale) F df(x, df1, df2) pf(q, df1, df2) qf(p, df1, df2) rf(n, df1, df2) T dt(x, df) pt(q, df) qt(p, df) rt(n, df) Chi-squared dchisq(x, df) pchi(q, df) qchisq(p, df) rchisq(n, df) 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. 6.3 Các hàm phân phối xác suất (probability distribution function) 6.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 4
  56. đượ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à: nk nk− Pknp()|,=− Cpk () 1 p , trong đó k == 0, 1, 2, . . . , n. Để hiểu định lí đó rõ ràng hơn, chúng ta sẽ xem qua qua vài ví dụ sau đây. Ví dụ 1: Hàm mật độ nhị phân (Binomial density probability function). Trong ví dụ trên, lớp học có 10 người, trong đó có 6 nữa. Nếu 3 bạn được chọn một cách ngẫu nhiên, xác suất mà chúng ta có 2 bạn nữ là bao nhiêu? Chúng ta có thể trả lời câu hỏi này một cách tương đối thủ công bằng cách xem xét tất cả các trường hợp có thể xảy ra. Mỗi lần chọn có 2 khả khăng (nam hay nữ), và 3 lần chọn, chúng ta có 23 = 8 trường hợp như sau. Bạn 1 Bạn 2 Bạn 3 Xác suất Nam Nam Nam (0.4)(0.4)(0.4) = 0.064 Nam Nam Nữ (0.4)(0.4)(0.6) = 0.096 Nam Nữ Nam (0.4)(0.6)(0.4) = 0.096 Nam Nữ Nữ (0.4)(0.6)(0.6) = 0.144 Nữ Nam Nam (0.6)(0.4)(0.4) = 0.096 Nữ Nam Nữ (0.6)(0.4)(0.6) = 0.144 Nữ Nữ Nam (0.6)(0.6)(0.4) = 0.144 Nữ Nữ Nữ (0.6)(0.6)(0.6) = 0.216 Tất cả các trường hợp 1.000 Chứng ta biết trước rằng trong nhóm 10 học sinh có 6 nữ, và do đó, xác suất nữ là 0.60. (Nói cách khác, xác suất chọn một bạn nam là 0.4). Do đó, xác suất mà tất cả 3 bạn được chọn đều là nam giới là: 0.4 x 0.4 x 0.4 = 0.064. Trong bảng trên, chúng ta thấy có 3 trường hợp mà trong đó có 2 bạn gái: đó là trường hợp Nam-Nữ-Nữ, Nữ-Nữ-Nam, và Nữ-Nam-Nữ, cả 3 đều có xác suất 0.144. Thành ra, xác suất chọn đúng 2 bạn nữ trong số 3 bạn được chọn là 3x0.144= 0.432. Trong R, có hàm dbinom(k, n, p) có thể giúp chúng ta tính công thức nk nk− Pknp()|,=− Cpk () 1 p một cách nhanh chóng. Trong trườ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 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à: 5
  57. > 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 table(b) b 0 1 2 3 4 5 6 7 8 9 10 6 45 147 192 229 169 105 68 23 13 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") Number of hypertensive patients Frequency 0 50 100 150 200 0246810 b 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 6
  58. 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). Ví dụ 4: Ứng dụng hàm phân phối nhị phân: Hai mươi khách hàng được mời uống hai loại bia A và B, và được hỏi họ thích bia nào. Kết quả cho thấy 16 người thích bia A. Vấn đề đặt ra là kết quả này có đủ để kết luận rằng bia A được nhiều người thích hơn bia B, hay là kết quả chỉ là do các yếu tố ngẫu nhiên gây nên? Chúng ta bắt đầu giải quyết vấn đề bằng cách giả thiết rằng nếu không có khác nhau, thì xác suất p=0.50 thích bia A và q=0.5 thích bia B. Nếu giả thiết này đúng, thì xác suất mà chúng ta quan sát 16 người trong số 20 người thích bia A là bao nhiêu. Chúng ta có thể tính xác suất này bằng R rất đơn giản: > 1- pbinom(15, 20, 0.5) [1] 0.005908966 Đáp số là xác suất 0.005 hay 0.5%. Nói cách khác, nếu quả thật hai bia giống nhau thì xác suất mà 16/20 người thích bia A chỉ 0.5%. Tức là, chúng ta có bằng chứng cho thấy khả năng bia A quả thật được nhiều người thích hơn bia B, chứ không phải do yếu tố ngẫu nhiên. Chú ý, chúng ta dùng 15 (thay vì 16), là bởi vì P(X ≥ 16) = 1 – P(X ≤ 15). Mà trong trường hợp ta đang bàn, P(X ≤ 15) = pbinom(15, 20, 0.5). 6.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ụ 5: 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ó 7
  59. 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 λ, : e−λ λ k PX()== k| λ k! e−221 Do đó, đáp số cho câu hỏi trên là: PX()===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) [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: PX()>=2 PX( =+ 3) PX( =+ 4) PX ( =+ 5) = 12− PX( ≤ ) = 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 6.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 8
  60. 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à:  2  2 1 ()x − µ PX()===− x|,µσ f() x exp 2  2σ 2πσ   Ví dụ 6: 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: Probability distribution of height in Vietnamese women f(height) 0.00 0.02 0.04 0.06 0.08 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 plot(height, dnorm(height, 156, 4.6), type="l", ylab=”f(height)”, xlab=”Height”, 9
  61. 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à: 2 1  ()160− 156  P(X = 160 | µ=156, σ=4.6) = exp − 2  4.6 2× 3.1416  24.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 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 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. 10
  62. Ví dụ 7: Ứ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ười trong quần thể này có áp suất máu bằng hoặc cao hơn 120 mmHg. 6.3.4 Hàm phân phối chuẩn chuẩn hóa (Standardized Normal distribution) 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 thường được viết tắt là: X ~ N(µ , σ2) Ở đâ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: 11
  63. Probability distribution of height in Vietnamese women f(z) 0.0 0.1 0.2 0.3 0.4 -4 -2 0 2 4 z Biểu đồ 3. Phân phối chuẩn hóa chiều cao ở phụ nữ Việt Nam. Biểu đồ trên được vẽ bằng hai lệnh sau đây: > height 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 pnorm(1.96) - pnorm(-1.96) [1] 0.9500042 12
  64. 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ụ 6 (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 qnorm(0.95, mean=0, sd=1) [1] 1.644854 Hay P(Z qnorm(0.975, mean=0, sd=1) [1] 1.959964 6.3.5 Hàm phân phối t, F và χ2 Các hàm phân phối t, F và χ2 trong thực tế là hàm của hàm phân phối chuẩn. Mối liên hệ và cách tính các hàm này có thể được mô tả bằng vài ghi chú sau đây: • Phân phối Chi bình phương (χ2). Phân phối χ2 xuất phát từ tổng bình n 2 phương của một biến phân phối chuẩn. Nếu nếu xi ~ N(0, 1), và gọi ux= ∑ i , thì i= u tuân theo luật phân phối Chi bình phương với bậc tự do n (thường viết tắt là 2 df). Nói theo ngôn ngữ toán, u ~ χn . 13
  65. Ví dụ 9: Tìm xác suất của một biến Chi bình phương, do đó, chỉ cần hai thông số u và n. Chẳng hạn như nếu chúng ta muốn tìm xác suất P(u=21, df=13), chỉ đơn giản dùng hàm pchisq như sau: > dchisq(21, 13) [1] 0.01977879 Tìm xác suất mà một biến số u nhỏ hơn 21 với bậc tự do 13 df. Tức là tìm P(u ≤ 21 | df=13) = ? > pchisq(21, 13) [1] 0.9270714 2 Cũng có thể nói kết quả trên cho biết P( χ13 qchisq(0.95, 15) [1] 24.99579 2 Nói cách khác, P( χ15 pchisq(21, 13, 5.4) [1] 0.6837649 2 Tức là, P( χ13,5.4 < 21) = 0.684. 14
  66. Tìm quantile của một trị số tương đương với 50% của một phân phối χ2 với 7 bậc tự do và thông số non-centrality bằng 3. > qchisq(0.5, 7, 3) [1] 9.180148 2 Do đó, P( χ7,3 1-pt(1.1, 6) [1] 0.1567481 Tức là, P(t6 > 1.1) = 1 – P(t6 qt(0.95, 15) [1] 1.753050 Nói cách khác, P(t19 1-pf(3.24, 3, 15, 5) [1] 0.3558721 Do đó, P(F3, 15, 5 > 3.24) = 1 - P(F3, 15,5 ≥ 3.24) = 0.355338. Với bậc tự do 3 và 15, tìm C sao cho P(F3, 15 > C) = 0.05. Lời giải của R là: 15
  67. > qf(1-0.05, 3, 15) [1] 3.287382 Nói cách khác, P(F3, 15 > 3.287382) = 1 – P(F3, 15 ≥ 3.287382) = 1 – 0.95 = 0.05 6.4 Mô phỏng (simulation) Trong phân tích thống kê, đôi khi vì hạn chế số mẫu chúng ta khó có thể ước tính một cách chính xác các thông số, và trong trường hợp bất định đó, chúng ta cần đến mô phỏng để biết được độ dao động của một hay nhiều thông số. Mô phỏng thường dựa vào các luật phân phối. Đây là một lĩnh vực khá phức tạp mà tôi không có ý định trình bày đầy đủ trong chương này. Ở đây, tôi chỉ trình bày một số mô hình mô phỏng mang tính minh họa để bạn đọc có thể dựa vào đó mà phát triển thêm. Ví dụ 11: Mô phỏng để chứng minh phương sai của số trung bình bằng phương sai chia cho n (var()X = σ 2 / n ). Chúng ta sẽ xem một biến số không liên tục với giá trị 1, 3 và 5 với xác suất như sau: x P(x) 1 0.60 3 0.30 5 0.10 Qua số liệu này, chúng ta biết rằng giá trị trung bình là (1x0.60)+(3x0.30)+(5x0.10) = 2.0 và phương sai (bạn đọc có thể tự tính) là 1.8. Bây giờ chúng ta sử dụng hai thông số này để thử mô phỏng 500 lần. Lệnh thứ nhất tạo ra 3 giá trị của x. Lệnh thứ hai nhập số xác suất cho từng giá trị của x. Lệnh sample yêu cầu R tạo nên 500 số ngẫu nhiên và cho vào đối tượng draws. x <- c(1, 3, 5) px <- c(0.6, 0.3, 0.1) draws <- sample(x, size=500, replace=T, prob=px) hist(draws, breaks=seq(1,5, by=0.25), main=”1000 draws”) 16
  68. 500 draws Frequency 0 50 100 150 200 250 300 12345 draws Từ luật phân phối xác suất chúng ta biết rằng tính trung bình sẽ có 60% lần có giá trị “1”, 30% có giá trị “2”, và 10% có giá trị “5”. Do đó, chúng ta kì vọng sẽ quan sát 300, 150 và 50 lần cho mỗi giá trị. Biểu đồ trên cho thấy phân phối các giá trị này gần với giá trị mà chúng ta kì vọng. Ngoài ra, chúng ta cũng biết rằng phương sai của biến số này là khoảng 1.8. Bây giờ chúng ta kiểm tra xem có đúng như kì vọng hay không: > var(draws) [1] 1.835671 Kết quả trên cho thấy phương sai của 500 mẫu là 1.836, tức không xa mấy so với giá trị kì vọng. Bây giờ chúng ta thử mô phỏng 500 giá trị trung bình x ( x là số trung bình của 4 số liệu mô phỏng) từ quần thể trên: > draws draws = matrix(draws, 4) > drawmeans = apply(draws, 2, mean) Lệnh thứ nhất và thứ hai tạo nên đối tượng tên là draws với 4 dòng, mỗi dòng có 500 giá trị từ luật phân phối trên. Nói cách khác, chúng ta có 4*500 = 2000 số. 500 số cũng có nghĩa là 500 cột: 1 đến 500. Tức mỗi cột có 4 số. Lệnh thứ ba tìm trị số trung bình cho mỗi cột. Lệnh này sẽ cho ra 500 số trung bình và chứa trong đối tượng drawmeans. Biểu đồ sau đây cho thấy phân phối của 500 số trung bình: > hist(drawmeans,breaks=seq(1,5,by=0.25), main=”1000 means of 4 draws”) 17
  69. 1000 means of 4 draws Frequency 0 50 100 150 12345 drawmeans Chúng ta thấy rằng phương sai của phân phối này nhỏ hơn. Thật ra, phương sai của 500 số trung bình này là 0.45. > var(drawmeans) [1] 0.4501112 Đây là giá trị tương đương với giá trị 0.45 mà chúng ta kì vọng từ công thức var()X ===σ 2 / 4 1.8/ 4 0.45. 6.4.1 Mô phỏng phân phối nhị phân Ví dụ 12: Mô phỏng mẫu từ một quần thể với luật phân phối nhị phân. Giả dụ chúng ta biết một quần thể có 20% người bị bệnh đái đường (xác suất p=0.2). Chúng ta muốn lấy mẫu từ quần thể này, mỗi mẫu có 20 đối tượng, và phương án chọn mẫu được lặp lại 100 lần: > bin bin [1] 4 4 5 3 2 2 3 2 5 4 3 6 7 3 4 4 1 5 3 5 3 4 4 5 1 4 4 4 4 3 2 4 2 2 5 4 5 [38] 7 3 5 3 3 4 3 2 4 5 2 4 5 5 4 2 2 2 8 5 5 5 3 4 5 7 4 3 6 4 6 6 8 8 3 3 1 [75] 1 4 4 2 3 9 7 4 4 0 0 8 6 9 3 1 4 5 6 4 5 3 2 4 3 2 Kết quả trên là số lần đầu, chúng ta sẽ có 4 người mắc bệnh; lần 2 cũng 4 người; lần 3 có 5 người mắc bệnh; v.v kết quả này có thể tóm lược trong một biểu đồ như sau: > hist(bin, xlab=”Number of diabetic patients”, ylab=”Number of samples”, main=”Distribution of the number of diabetic patients”) 18
  70. Distribution of the number of diabetic patients Number samples of 0 5 10 15 20 25 02468 Number of diabetic patients > mean(bin) [1] 3.97 Đúng như chúng ta kì vọng, vì chọn mỗi lần 20 đối tượng và xác suất 20%, nên chúng ta tiên đoán trung bình sẽ có 4 bệnh nhân đái đường. 6.4.2 Mô phỏng phân phối Poisson Ví dụ 13: Mô phỏng mẫu từ một quần thể với luật phân phối Poisson. Trong ví dụ sau đây, chúng ta mô phỏng 100 mẫu từ một quần thể tuân theo luật phân phối Poisson với trung bình λ=3: > pois pois > pois [1] 4 3 2 4 2 3 4 4 0 7 5 0 3 3 4 2 2 6 1 4 2 3 3 5 4 2 1 4 0 2 1 5 1 2 2 2 6 [38] 1 3 6 3 3 5 4 3 2 2 5 3 3 3 1 4 7 3 4 3 2 6 1 4 1 0 5 2 2 2 3 6 8 4 4 1 4 [75] 1 0 0 4 3 3 2 3 3 3 4 1 5 4 4 1 3 1 6 4 4 4 2 2 2 4 Và mật độ phân phối: 19
  71. Histogram of pois Frequency 0 5 10 15 20 02468 pois Phân phối Poisson và phân phối mũ. Trong ví dụ sau đây, chúng ta mô phỏng thời gian bệnh nhân đến một bệnh viện. Biết rằng bệnh nhân đến bệnh viện một cách ngẫu nhiên theo luật phân phối Poisson, với trung bình 15 bệnh nhân cho mỗi 150 phút. Có thể chứng minh dễ dàng rằng thời gian giữa hai bệnh nhân đến bệnh viện tuân theo luật phân phối mũ. Chúng ta muốn biết thời gian mà bệnh nhân ghé bệnh viện; do đó, chúng ta mô phỏng 15 thời gian giữa hai bệnh nhân từ luật phân phối mũ với tỉ lệ 15/150 = 0.1 mỗi phút. Các lệnh sau đây đáp ứng yêu cầu đó: # Tạo thời gian đến bệnh viện > appoint times times [1] 37 5 8 10 24 5 1 7 8 6 12 6 3 25 15 6.4.3 Mô phỏng phân phối χ2, t, F Cách mô phỏng trên đây còn có thể áp dụng cho các luật phân phối khác như nhị phân âm (negative binomial distribution với rnbinom), gamma (rgamma), beta (rbeta), Chi bình phương (rchisq), hàm mũ (rexp), t (rt), F (rf), v.v Các thông số cho các hàm mô phỏng này có thể tìm trong phần đầu của chương. Các lệnh sau đây sẽ minh họa các luật phân phối thông thường đó: • Phân phối Chi bình phương với một số bậc tự do: > curve(dchisq(x, 1), xlim=c(0,10), ylim=c(0,0.6), col="red", lwd=3) > curve(dchisq(x, 2), add=T, col="green", lwd=3) > curve(dchisq(x, 3), add=T, col="blue", lwd=3) > curve(dchisq(x, 5), add=T, col="orange", lwd=3) > abline(h=0, lty=3) 20
  72. > abline(v=0, lty=3) > legend(par("usr")[2], par("usr")[4], xjust=1, c("df=1", "df=2", "df=3", "df=5"), lwd=3, lty=1, col=c("red", "green", "blue", "orange")) df=1 df=2 df=3 df=5 dchisq(x, 1) dchisq(x, 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0246810 x Biểu đồ 4. Phân phối Chi bình phương với bậc tự do =1, 2, 3, 5. • Phân phối t: > curve(dt(x, 1), xlim=c(-3,3), ylim=c(0,0.4), col="red", lwd=3) > curve(dt(x, 2), add=T, col="blue", lwd=3) > curve(dt(x, 5), add=T, col="green", lwd=3) > curve(dt(x, 10), add=T, col="orange", lwd=3) > curve(dnorm(x), add=T, lwd=4, lty=3) > title(main=“Student T distributions”) > legend(par("usr")[2], par("usr")[4], xjust=1, c("df=1", "df=2", "df=5", "df=10", “Normal distribution”), lwd=c(2,2,2,2,2), lty=c(1,1,1,1,3), col=c("red", "blue", "green", "orange", par(“fg”))) 21
  73. Student T distributions df=1 df=2 df=5 df=10 Normal distribution dt(x, 1) dt(x, 0.0 0.1 0.2 0.3 0.4 -3 -2 -1 0 1 2 3 x Biểu đồ 5. Phân phối t với bậc tự do =1, 2, 5, 10 so sánh với phân phối chuẩn. • Phân phối F: > curve(df(x,1,1), xlim=c(0,2), ylim=c(0,0.8), lwd=3) > curve(df(x,3,1), add=T) > curve(df(x,6,1), add=T, lwd=3) > curve(df(x,3,3), add=T, col="red") > curve(df(x,6,3), add=T, col="red", lwd=3) > curve(df(x,3,6), add=T, col="blue") > curve(df(x,6,6), add=T, col="blue", lwd=3) > title(main=“Fisher F distributions”) > legend(par("usr")[2], par("usr")[4], xjust=1, c("df=1,1", "df=3,1", "df=6,1", "df=3,3", “df=6,3”, “df=3,6”, df=”6,6”), lwd=c(1,1,3,1,3,1,3), lty=c(2,1,1,1,1,1,1), col=c(par("fg"), par("fg"), par("fg"), “red”, “blue”, “blue”)) 22
  74. Fisher F distributions df=1,1 df=3,1 df=6,1 df=3,3 df=6,3 df=3,6 6,6 df(x, 1, 1) 1, df(x, 0.00.20.40.60.8 0.00.51.01.52.0 x Biểu đồ 6. Phân phối F với nhiều bậc tự do khác nhau. • Phân phối gamma: > curve( dgamma(x,1,1), xlim=c(0,5) ) > curve( dgamma(x,2,1), add=T, col='red' ) > curve( dgamma(x,3,1), add=T, col='green' ) > curve( dgamma(x,4,1), add=T, col='blue' ) > curve( dgamma(x,5,1), add=T, col='orange' ) > title(main="Gamma probability distribution function") > legend(par('usr')[2], par('usr')[4], xjust=1, c('k=1 (Exponential distribution)', 'k=2', 'k=3', 'k=4', 'k=5'), lwd=1, lty=1, col=c(par('fg'), 'red', 'green', 'blue', 'orange') ) Gamma probability distribution function k=1 (Exponential distribution) k=2 k=3 k=4 k=5 dgamma(x, 1, 1) 1, dgamma(x, 0.0 0.2 0.4 0.6 0.8 1.0 012345 x 23
  75. Biểu đồ 7. Phân phối Gamma với nhiều hình dạng. • Phân phối beta: > curve( dbeta(x,1,1), xlim=c(0,1), ylim=c(0,4) ) > curve( dbeta(x,2,1), add=T, col='red' ) > curve( dbeta(x,3,1), add=T, col='green' ) > curve( dbeta(x,4,1), add=T, col='blue' ) > curve( dbeta(x,2,2), add=T, lty=2, lwd=2, col='red' ) > curve( dbeta(x,3,2), add=T, lty=2, lwd=2, col='green' ) > curve( dbeta(x,4,2), add=T, lty=2, lwd=2, col='blue' ) > curve( dbeta(x,2,3), add=T, lty=3, lwd=3, col='red' ) > curve( dbeta(x,3,3), add=T, lty=3, lwd=3, col='green' ) > curve( dbeta(x,4,3), add=T, lty=3, lwd=3, col='blue' ) > title(main="Beta distribution") > legend(par('usr')[1], par('usr')[4], xjust=0, c('(1,1)', '(2,1)', '(3,1)', '(4,1)', '(2,2)', '(3,2)', '(4,2)', '(2,3)', '(3,3)', '(4,3)' ), lwd=1, #c(1,1,1,1, 2,2,2, 3,3,3), lty=c(1,1,1,1, 2,2,2, 3,3,3), col=c(par('fg'), 'red', 'green', 'blue', 'red', 'green', 'blue', 'red', 'green', 'blue' )) Beta distribution (1,1) (2,1) (3,1) (4,1) (2,2) (3,2) (4,2) (2,3) (3,3) (4,3) dbeta(x, 1, 1) 1, dbeta(x, 01234 0.0 0.2 0.4 0.6 0.8 1.0 x Biểu đồ 8. Phân phối beta với nhiều hình dạng. • Phân phối Weibull: > curve(dexp(x), xlim=c(0,3), ylim=c(0,2)) > curve(dweibull(x,1), lty=3, lwd=3, add=T) > curve(dweibull(x,2), col='red', add=T) 24
  76. > curve(dweibull(x,.8), col='blue', add=T) > title(main="Weibull Probability Distribution Function") > legend(par('usr')[2], par('usr')[4], xjust=1, c('Exponential', 'Weibull, shape=1', 'Weibull, shape=2', 'Weibull, shape=.8'), lwd=c(1,3,1,1), lty=c(1,3,1,1), col=c(par("fg"), par("fg"), 'red', 'blue')) Weibull Probability Distribution Function Exponential Weibull, shape=1 Weibull, shape=2 Weibull, shape=.8 dexp(x) 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 x Biểu đồ 9. Phân phối Weibull. • Phân phối Cauchy: > curve(dcauchy(x),xlim=c(-5,5), ylim=c(0,.5), lwd=3) > curve(dnorm(x), add=T, col='red', lty=2) > legend(par('usr')[2], par('usr')[4], xjust=1, c('Cauchy distribution', 'Gaussian distribution'), lwd=c(3,1), lty=c(1,2), col=c(par("fg"), 'red')) 25
  77. Cauchy distribution Gaussian distribution dcauchy(x) 0.0 0.1 0.2 0.3 0.4 0.5 -4 -2 0 2 4 x Biểu đồ 9. Phân phối Cauchy so sánh với phân phối chuẩn. 6.5 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ụ: 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 26
  78. 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) [1] 31 44 6 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. 27