1. Thống kê mô tả(descriptive statistíc, summary) Để minh họa cho việc áp dung 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à mydata. Trong nghiên cứu này , ngoài các chỉ số liên quan đến category, city, Country, v.v tôi còn đo lường các liên quan đến Profit.Ratio, Profit.per.Order, Sales.Forecast. Có 62 chỉ số nghiên cứu. Dữ liệu này chứa trong directory: E:\quyennv\RStudio. 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: file .csv đính đính kèm ở dưới:
# xem xét các cột trong dữ liệu > names(mydata) > mydata
Cho một biến x1, x2, x3, ..., xn chúng ta có thể tính toán một chỉ số thông kê mô tả như sau:
Để tìm giá trị trung bình của Quantity, chúng ta chỉ dùng lệnh: > mean(Quantity) [1] 3.74
Hay phương sai và độ lệnh chuẩn của Quantity: > var(Quantity) [1] 4.234747
> sd(Quantity) [1] 2.05785
Tuy nhiên , R có lệnh Summary, có thể cho chúng ta biết tất cả thông tin thống kê về một biến số: > summary(Quantity) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 2.00 3.00 3.74 5.00 9.00 Trong đó, "1st Qu, 3rd Qu" có nghĩa là first quartile (tương đương với giá trị 25%) và third quartile (tương đương với giá trị 75%) của một biến số. First quartile = 2 coa nghĩa là 25% đối tượng nghiên cứu có Quantity(số lượng) bằng hoặc nhỏ hơn 2. Tương tự, third quartile = 5 có nghĩa là 75% có Quantity bằng hoặc thấp hơn 5. Tất nhiên số trung vị (median) 3.74 cũng có nghĩa là 50%.
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 bất cứ khi nào ta muốn > desc(Quantity) MEAN SD SE 3.740000 2.057850 0.205785
Để xem phân phối của hormones và chỉ số cùng một lúc, chúng ta có thể vẽ đồ thị cho biến số. Trước hết, chia màng hình thành 4 cửa sổ(với 2 dòng và 2 cột); Sau đó lần lượt vẽ: > op <- par(mfrow=c(2,2)) > hist(Quantity) > hist(Days.to.Ship.Actual) > hist(Days.to.Ship.Scheduled)
2. 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.
2.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ế. Chúng ta có dữ liệu gồm 8 cột (hay 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:
2.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: > barplot(table(sex), main=”Frequency of males and females”)
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 > barplot(sex.freq, horiz = TRUE, col = rainbow(length(sex.freq)), main=”Frequency of males and females”)
2.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] (67.3,80] 19 24 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:
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”)
2.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)) Ví dụ pie(table(cut(age,5)))
2.5 Biểu đồ cho một biến số liên tục: stripchart và hist 2.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”)
2.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")
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 > plot(density(age),add=TRUE)
2.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:
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")
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: > boxplot(tc~sex, horizontal=TRUE, main="Box plot of total cholesterol", ylab="mg/L", col = "pink")
2.7 Phân tích biểu đồ cho hai biến liên tục 2.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)
Chúng ta mu 8 ố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)) > plot(hdl, tc, pch=ifelse(sex=="Nam", “M”, “F”))
1. Hoán Vị(permutation) Chúng ta biết 3! = 3.2.1 = 6 và 0!=1. Nói chung, công thức tính số hoán vị cho 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ổ hợp tính bằng hàm choose(n,k) Ví dụ choose(5,2) = 10
3. Biến ngẩu nhiên và biến 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 có thể có. Các hàm phân phối (distribution function) là hàm mô tả các biến đó một cách hệ thống. “Có hệ thống” ở đây có nghĩa là theo một 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, chúng ta sẽ em 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 đinh sẵ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), R NDH 19 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:
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.
a. Hàm phân phối nhị phân Như tên gọi, hàm phân phối nhị phân có 2 giá trị: Nam/ Nữ, có/ không,v.v... Hàm nhị phân được phát biểu bằng định lí sau: Nếu một thử nghiêm được thực hành n lần, mổi lần cho ra kết quả(thành công hoặc thất bại), và 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à: P(k | n, p) = C(n,k)p^k(1-p)^n-k, trong đó k=0, 1,2,3,..., n. Trong R ta có hàm dbinom(k, n, p) có thể giúp ta tính nhanh một cách nhanh chóng. Ví dụ: > dbinom(2, 3, 0.60) [1] 0.432
b. 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!
Do đó, đáp số cho câu hỏi trên là: e ^-1 /2! = 0,1839 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 > dpois(2,1)
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 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
c. 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 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ó dạng:
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
Hàm xác suất chuẩn tích lũy (cumulative normal probability function). Vì hiề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(150 ≤ X ≤ 160), 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:
d. Hàm phân phối chuẩn chuẩn hóa (Standardized Normal distribution) 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 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=). Ví dụ: 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
