0.0008. Trong chương này, chúng ta sẽ làm quen với một số lệnh trong R để tiến hành
những tính toán đơn giản trên.
9.1 Thống kê mô tả (descriptive statistics, summary)
Để minh họa cho việc áp dụng R vào thống kê mô tả, tôi sẽ sử dụng một dữ liệu
nghiên cứu có tên là igfdata. Trong nghiên cứu này, ngoài các chỉ số liên quan đến
giới tính, độ tuổi, trọng lượng và chiều cao, chúng tôi đo lường các hormone liên quan
đến tình trạng tăng trưởng như igfi, igfbp3, als, và các markers liên quan đến
sự chuyển hóa của xương pinp, ictp và pinp. Có 100 đối tượng nghiên cứu. Dữ
liệu này được chứa trong directory c:\works\stats. Trước hết, chúng ta cần phải
nhập dữ liệu vào R với những lệnh sau đây (các câu chữ theo sau dấu # là những chú
thích để bạn đọc theo dõi):
> options(width=100)
# chuyển directory
> setwd("c:/works/stats")
# đọc dữ liệu vào R
> igfdata <- read.table("igf.txt", header=TRUE, na.strings=".")
> attach(igfdata)
# xem xét các cột số trong dữ liệu
> names(igfdata)
[1] "id" "sex" "age" "weight" "height" "ethnicity"
[7] "igfi" "igfbp3" "als" "pinp" "ictp" "p3np"
> igfdata
id sex age weight height ethnicity igfi igfbp3 als pinp ictp p3np
1 1 Female 15 42 162 Asian 189.000 4.00000 323.667 353.970 11.2867 8.3367
2 2 Male 16 44 160 Caucasian 160.000 3.75000 333.750 375.885 10.4300 6.7450
3 3 Female 15 43 157 Asian 146.833 3.43333 248.333 199.507 8.3633 12.5000
4 4 Female 15 42 155 Asian 185.500 3.40000 251.000 483.607 13.3300 14.2767
5 5 Female 16 47 167 Asian 192.333 4.23333 322.000 105.430 7.9233 4.5033
6 6 Female 25 45 160 Asian 110.000 3.50000 284.667 76.487 4.9833 4.9367
7 7 Female 19 45 161 Asian 157.000 3.20000 274.000 75.880 6.3500 5.3200
8 8 Female 18 43 153 Asian 146.000 3.40000 303.000 86.360 7.3700 4.6700
9 9 Female 15 41 149 Asian 197.667 3.56667 308.500 254.803 11.8700 6.8200
10 10 Female 24 45 157 African 148.000 3.40000 273.000 44.720 3.7400 6.1600
97 97 Female 17 54 168 Caucasian 204.667 4.96667 441.333 64.130 5.1600 4.4367
98 98 Male 18 55 169 Asian 178.667 3.86667 273.000 185.913 7.5267 8.8333
99 99 Female 18 48 151 Asian 237.000 3.46667 324.333 105.127 5.9867 5.6600
100 100 Male 15 54 168 Asian 130.000 2.70000 259.333 325.840 10.2767 6.5933
Trên đây chỉ là một phần số liệu trong số 100 đối tượng.
Cho một biến số
123
, , , ,
n
x xx x chúng ta có thể tính toán một số chỉ số thống kê mô tả
như sau:
Lí thuyết
Hàm R
Số trung bình:
x
n
x
i
i
n
=
=
∑
1
1
.
mean(x)
Phương sai:
()
∑
−
−
=
=
n
i
i
xx
n
s
1
2
2
1
1
var(x)
Độ lệch chuẩn:
2
ss=
sd(x)
Sai số chuẩn (standard error):
s
SE
n
=
Không có
Trị số thấp nhất
min(x)
Trị số cao nhất
max(x)
Toàn cự (range)
range(x)
Ví dụ 1
: Để tìm giá trị trung bình của độ tuổi, chúng ta chỉ đơn giản lệnh:
> mean(age)
[1] 19.17
Hay phương sai và độc lệch chuẩn của tuổi:
> var(age)
[1] 15.33444
> sd(age)
[1] 3.915922
Tuy nhiên,
R
có lệnh
summary
có thể cho chúng ta tất cả thông tin thống kê về một biến
số:
> summary(age)
Min. 1st Qu. Median Mean 3rd Qu. Max.
13.00 16.00 19.00 19.17 21.25 34.00
Nói chung, kết quả này đơn giản và các viết tắt cũng có thể dễ hiểu. Chú ý, trong
kết quả trên, có hai chỉ số “
1st Qu
” và “
3rd Qu
” có nghĩa là first quartile (tương
đương với vị trí 25%) và third quartile (tương đương với vị trí 75%) của một biến số.
First quartile = 16 có nghĩa là 25% đối tượng nghiên cứu có độ tuổi bằng hoặc nhỏ hơn
16 tuổi. Tương tự, Third quartile = 34 có nghĩa là 75% đối tượng có độ tuổi bằng hoặc
thấp hơn 34 tuổi. Tất nhiên số trung vị (median) 19 cũng có nghĩa là 50% đối tượng có
độ tuổi 19 trở xuống (hay 19 tuổi trở lên).
R
không có hàm tính sai số chuẩn, và trong hàm summary,
R
cũng không cung
cấp độ lệch chuẩn. Để có các số này, chúng ta có thể tự viết một hàm đơn giản (hãy gọi
là
desc
) như sau:
desc <- function(x)
{
av <- mean(x)
sd <- sd(x)
se <- sd/sqrt(length(x))
c(MEAN=av, SD=sd, SE=se)
}
Và có thể gọi hàm này để tính bất cứ biến nào chúng ta muốn, như tính biến
als
sau
đây:
> desc(als)
MEAN SD SE
301.841120 58.987189 5.898719
Để có một “quang cảnh” chung về dữ liệu
igfdata
chúng ta chỉ đơn giản lệnh
summary
như sau:
> summary(igfdata)
id sex age weight height ethnicity
Min. : 1.00 Female:69 Min. :13.00 Min. :41.00 Min. :149.0 African : 8
1st Qu.: 25.75 Male :31 1st Qu.:16.00 1st Qu.:47.00 1st Qu.:157.0 Asian :60
Median : 50.50 Median :19.00 Median :50.00 Median :162.0 Caucasian:30
Mean : 50.50 Mean :19.17 Mean :49.91 Mean :163.1 Others : 2
3rd Qu.: 75.25 3rd Qu.:21.25 3rd Qu.:53.00 3rd Qu.:168.0
Max. :100.00 Max. :34.00 Max. :60.00 Max. :196.0
igfi igfbp3 als pinp ictp
Min. : 85.71 Min. :2.000 Min. :192.7 Min. : 26.74 Min. : 2.697
1st Qu.:137.17 1st Qu.:3.292 1st Qu.:256.8 1st Qu.: 68.10 1st Qu.: 4.878
Median :161.50 Median :3.550 Median :292.5 Median :103.26 Median : 6.338
Mean :165.59 Mean :3.617 Mean :301.8 Mean :167.17 Mean : 7.420
3rd Qu.:186.46 3rd Qu.:3.875 3rd Qu.:331.2 3rd Qu.:196.45 3rd Qu.: 8.423
Max. :427.00 Max. :5.233 Max. :471.7 Max. :742.68 Max. :21.237
p3np
Min. : 2.343
1st Qu.: 4.433
Median : 5.445
Mean : 6.341
3rd Qu.: 7.150
Max. :16.303
R
tính toán tất cả các biến số nào có thể tính toán được! Thành ra, ngay cả cột
id
(tức mã số của đối tượng nghiên cứu)
R
cũng tính luôn! (và chúng ta biết kết quả của cột
id
chẳng có ý nghĩa thống kê gì). Đối với các biến số mang tính phân loại như
sex
và
ethnicity
(sắc tộc) thì
R
chỉ báo cáo tần số cho mỗi nhóm.
Kết quả trên cho tất cả đối tượng nghiên cứu. Nếu chúng ta muốn kết quả cho
từng nhóm nam và nữ riêng biệt, hàm
by
trong
R
rất hữu dụng. Trong lệnh sau đây,
chúng ta yêu cầu
R
tóm lược dữ liệu
igfdata
theo
sex
.
> by(igfdata, sex, summary)
sex: Female
id sex age weight height
Min. : 1.0 Female:69 Min. :13.00 Min. :41.00 Min. :149.0
1st Qu.:21.0 Male : 0 1st Qu.:17.00 1st Qu.:47.00 1st Qu.:156.0
Median :47.0 Median :19.00 Median :50.00 Median :162.0
Mean :48.2 Mean :19.59 Mean :49.35 Mean :161.9
3rd Qu.:75.0 3rd Qu.:22.00 3rd Qu.:52.00 3rd Qu.:166.0
Max. :99.0 Max. :34.00 Max. :60.00 Max. :196.0
ethnicity igfi igfbp3 als
African : 4 Min. : 85.71 Min. :2.767 Min. :204.3
Asian :43 1st Qu.:136.67 1st Qu.:3.333 1st Qu.:263.8
Caucasian:22 Median :163.33 Median :3.567 Median :302.7
Others : 0 Mean :167.97 Mean :3.695 Mean :311.5
3rd Qu.:186.17 3rd Qu.:3.933 3rd Qu.:361.7
Max. :427.00 Max. :5.233 Max. :471.7
pinp ictp p3np
Min. : 26.74 Min. : 2.697 Min. : 2.343
1st Qu.: 62.75 1st Qu.: 4.717 1st Qu.: 4.337
Median : 78.50 Median : 5.537 Median : 5.143
Mean :108.74 Mean : 6.183 Mean : 5.643
3rd Qu.:115.26 3rd Qu.: 7.320 3rd Qu.: 6.143
Max. :502.05 Max. :13.633 Max. :14.420
sex: Male
id sex age weight height
Min. : 2.00 Female: 0 Min. :14.00 Min. :44.00 Min. :155.0
1st Qu.: 34.50 Male :31 1st Qu.:15.00 1st Qu.:48.50 1st Qu.:161.5
Median : 56.00 Median :17.00 Median :51.00 Median :164.0
Mean : 55.61 Mean :18.23 Mean :51.16 Mean :165.6
3rd Qu.: 75.00 3rd Qu.:20.00 3rd Qu.:53.50 3rd Qu.:169.0
Max. :100.00 Max. :27.00 Max. :59.00 Max. :191.0
ethnicity igfi igfbp3 als
African : 4 Min. : 94.67 Min. :2.000 Min. :192.7
Asian :17 1st Qu.:138.67 1st Qu.:3.183 1st Qu.:249.8
Caucasian: 8 Median :160.00 Median :3.500 Median :276.0
Others : 2 Mean :160.29 Mean :3.443 Mean :280.2
3rd Qu.:183.00 3rd Qu.:3.775 3rd Qu.:311.3
Max. :274.00 Max. :4.500 Max. :388.7
pinp ictp p3np
Min. : 56.28 Min. : 3.650 Min. : 3.390
1st Qu.:135.07 1st Qu.: 6.900 1st Qu.: 5.375
Median :245.92 Median : 9.513 Median : 7.140
Mean :297.21 Mean :10.173 Mean : 7.895
3rd Qu.:450.38 3rd Qu.:13.517 3rd Qu.:10.010
Max. :742.68 Max. :21.237 Max. :16.303
Để xem qua phân phối của các hormones và chỉ số sinh hóa cùng một lúc, chúng
ta có thể vẽ đồ thị cho tất cả 6 biến số. Trước hết, chia màn ảnh thành 6 cửa sổ (với 2
dòng và 3 cột); sau đó lần lượt vẽ:
> op <- par(mfrow=c(2,3))
> hist(igfi)
> hist(igfbp3)
> hist(als)
> hist(pinp)
> hist(ictp)
> hist(p3np)
Histogram of igfi
igf i
Frequency
100 200 300 400
0 10203040
Histogram of igfbp3
igf bp3
Frequency
2.0 3.0 4.0 5.0
0 10203040
Histogram of a ls
als
Frequency
150 250 350 450
0 102030
Histogram of pinp
pinp
Frequency
0 200 400 600 800
01020304050
Histogram of ictp
ic tp
Frequency
5101520
0102030
Histogram of p3np
p3np
Frequency
51015
0 10203040
9.2 Kiểm định xem một biến có phải phân phối chuẩn
Trong phân tích thống kê, phần lớn các phép tính dựa vào giả định biến số phải là
một biến số phân phối chuẩn (normal distribution). Do đó, một trong những việc quan
trọng khi xem xét dữ kiện là phải kiểm định giả thiết phân phối chuẩn của một biến số.
Trong đồ thị trên, chúng ta thấy các biến số như
igfi, pinp, ictp
và
p3np
có vẻ
tập trung vào các giá trị thấp và không cân đối, tức dấu hiệu của một sự phân phối không
chuẩn.
Để kiểm định nghiêm chỉnh, chúng ta cần phải sử dụng kiểm định thống kê có tên
là “Shapiro test” và trong R gọi là hàm
shapiro.test
. Chẳng hạn như kiểm định giả
thiết phân phối chuẩn của biến số
pinp
,
> shapiro.test(pinp)
Shapiro-Wilk normality test
data: pinp
W = 0.748, p-value = 8.314e-12
Vì trị số p (p-value) thấp hơn 0.05, chúng ta có thể kết luận rằng biến số
pinp
không đáp
ứng luật phân phối chuẩn.
Nhưng với biến số
weight
(trọng lương cơ thể) thì kiểm định này cho biết đây là một
biến số tuân theo luật phân phối chuẩn vì trị số p > 0.05.
> shapiro.test(weight)
Shapiro-Wilk normality test
data: weight
W = 0.9887, p-value = 0.5587
Thật ra, kết quả trên cũng phù hợp với đồ thị của
weight
:
> hist(weight)
Histogram of weight
weight
Frequency
40 45 50 55 60
051015
9.3 Thống kê mô tả theo từng nhóm
Nếu chúng ta muốn tính trung bình của một biến số như igfi cho mỗi nhóm nam
và nữ giới, hàm
tapply
trong
R
có thể dùng cho việc này:
> tapply(igfi, list(sex), mean)
Female Male
167.9741 160.2903
Trong lệnh trên,
igfi
là biến số chúng ta cần tính, biến số phân nhóm là
sex
, và chỉ số
thống kê chúng ta muốn là trung bình (
mean
). Qua kết quả trên, chúng ta thấy số trung
bình của
igfi
cho nữ giới (167.97) cao hơn nam giới (160.29).
Nhưng nếu chúng ta muốn tính cho từng giới tính và sắc tộc, chúng ta chỉ cần thêm một
biến số trong hàm
list
:
> tapply(igfi, list(ethnicity, sex), mean)
Female Male
African 145.1252 120.9168
Asian 165.6589 160.4999
Caucasian 176.6536 169.4790
Others NA 200.5000
Trong kết quả trên,
NA
có nghĩa là “not available”, tức không có số liệu cho phụ nữ trong
các sắc tộc “others”.
9.4 Kiểm định t (t.test)
Kiểm định t dựa vào giả thiết phân phối chuẩn. Có hai loại kiểm định t: kiểm
định t cho một mẫu (one-sample t-test), và kiểm định t cho hai mẫu (two-sample t-test).
Kiểm định t một mẫu nằm trả lời câu hỏi dữ liệu từ một mẫu có phải thật sự bằng một
thông số nào đó hay không. Còn kiểm định t hai mẫu thì nhằm trả lời câu hỏi hai mẫu có
cùng một luật phân phối, hay cụ thể hơn là hai mẫu có thật sự có cùng trị số trung bình
hay không. Tôi sẽ lần lượt minh họa hai kiểm định này qua số liệu
igfdata
trên.
9.1.1 Kiểm định t một mẫu
Ví dụ 2
. Qua phân tích trên, chúng ta thấy tuổi trung bình của 100 đối tượng
trong nghiên cứu này là 19.17 tuổi. Chẳng hạn như trong quần thể này, trước đây chúng
ta biết rằng tuổi trung bình là 30 tuổi. Vấn đề đặt ra là có phải mẫu mà chúng ta có được
có đại diện cho quần thể hay không. Nói cách khác, chúng ta muốn biết giá trị trung bình
19.17 có thật sự khác với giá trị trung bình 30 hay không.
Để trả lời câu hỏi này, chúng ta sử dụng kiểm định t. Theo lí thuyết thống kê,
kiểm định t được định nghĩa bằng công thức sau đây:
/
x
t
sn
µ
−
=
Trong đó,
x
là giá trị trung bình của mẫu,
µ
là trung bình theo giả thiết (trong trường
hợp này, 30),
s
là độ lệch chuẩn, và
n
là số lượng mẫu (100). Nếu giá trị t cao hơn giá trị
lí thuyết theo phân phối t ở một tiêu chuẩn có ý nghĩa như 5% chẳng hạn thì chúng ta có
lí do để phát biểu khác biệt có ý nghĩa thống kê. Giá trị này cho mẫu 100 có thể tính toán
bằng hàm
qt
của
R
như sau:
> qt(0.95, 100)
[1] 1.660234
Nhưng có một cách tính toán nhanh gọn hơn để trả lời câu hỏi trên, bằng cách dùng hàm
t.test
như sau:
> t.test(age, mu=30)
One Sample t-test
data: age
t = -27.6563, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 30
95 percent confidence interval:
18.39300 19.94700
sample estimates:
mean of x
19.17
Trong lệnh trên
age
là biến số chúng ta cần kiểm định, và
mu=30
là giá trị giả thiết.
R
trình bày trị số
t = -27.66
, với 99 bậc tự do, và trị số p < 2.2e-16 (tức rất thấp).
R
cũng cho biết độ tin cậy 95% của age là từ 18.4 tuổi đến 19.9 tuổi (30 tuổi nằm quá ngoài
khoảng tin cậy này). Nói cách khác, chúng ta có lí do để phát biểu rằng độ tuổi trung
bình trong mẫu này thật sự thấp hơn độ tuổi trung bình của quần thể.
9.4.2 Kiểm định t hai mẫu
Ví dụ 3
. Qua phân tích mô tả trên (phầm
summary
) chúng ta thấy phụ nữ có độ
hormone
igfi
cao hơn nam giới (167.97 và 160.29). Câu hỏi đặt ra là có phải thật sự đó
là một khác biệt có hệ thống hay do các yếu tố ngẫu nhiên gây nên. Trả lời câu hỏi này,
chúng ta cần xem xét mức độ khác biệt trung bình giữa hai nhóm và độ lệch chuẩn của độ
khác biệt.
21
x x
t
SED
−
=
Trong đó
1
x
và
2
x
là số trung bình của hai nhóm nam và nữ, và
SED
là độ lệch chuẩn
của (
1
x
-
2
x
) . Thực ra,
SED
có thể ước tính bằng công thức:
22
12
SED SE SE=+
Trong đó
1
SE
và
2
SE
là sai số chuẩn (standard error) của hai nhóm nam và nữ. Theo lí
thuyết xác suất,
t
tuân theo luật phân phối
t
với bậc tự do
12
2
nn+ −
, trong đó
n
1
và
n
2
là
số mẫu của hai nhóm. Chúng ta có thể dùng
R
để trả lời câu hỏi trên bằng hàm
t.test
như sau:
> t.test(igfi~ sex)
Welch Two Sample t-test
data: igfi by sex
t = 0.8412, df = 88.329, p-value = 0.4025
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10.46855 25.83627
sample estimates:
mean in group Female mean in group Male
167.9741 160.2903
R
trình bày các giá trị quan trọng trước hết:
t = 0.8412, df = 88.329, p-value = 0.4025
df
là bậc tự do. Trị số p = 0.4025 cho thấy mức độ khác biệt giữa hai nhóm nam và nữ
không có ý nghĩa thống kê (vì cao hơn 0.05 hay 5%).
95 percent confidence interval:
-10.46855 25.83627
là khoảng tin cậy 95% về độ khác biệt giữa hai nhóm. Kết quả tính toán trên cho biết độ
igf
ở nữ giới có thể thấp hơn nam giới 10.5 ng/L hoặc cao hơn nam giới khoảng 25.8
ng/L. Vì độ khác biệt quá lớn và đó là thêm bằng chứng cho thấy không có khác biệt có
ý nghĩa thống kê giữa hai nhóm.
Kiểm định trên dựa vào giả thiết hai nhóm nam và nữ có khác phương sai. Nếu
chúng ta có lí do đề cho rằng hai nhóm có cùng phương sai, chúng ta chỉ thay đổi một
thông số trong hàm t với
var.equal=TRUE
như sau:
> t.test(igfi~ sex, var.equal=TRUE)
Two Sample t-test
data: igfi by sex
t = 0.7071, df = 98, p-value = 0.4812
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-13.88137 29.24909
sample estimates:
mean in group Female mean in group Male
167.9741 160.2903
Về mặc số, kết quả phân tích trên có khác chút ít so với kết quả phân tích dựa vào giả
định hai phương sai khác nhau, nhưng trị số p cũng đi đến một kết luận rằng độ khác biệt
giữa hai nhóm không có ý nghĩa thống kê.
9.5 So sánh phương sai (var.test)
Bây giờ chúng ta thử kiểm định xem phương sai giữa hai nhóm có khác nhau không. Để
tiến hành phân tích, chúng ta chỉ cần lệnh:
> var.test(igfi ~ sex)
F test to compare two variances
data: igfi by sex
F = 2.6274, num df = 68, denom df = 30, p-value = 0.004529
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.366187 4.691336
sample estimates:
ratio of variances
2.627396
Kết quả trên cho thấy độ khác biệt về phương sai giữa hai nhóm cao 2.62 lần. Trị số p =
0.0045 cho thấy phương sai giữa hai nhóm khác nhau có ý nghĩa thống kê. Như vậy,
chúng ta chấp nhận kết quả phân tích của hàm
t.test(igfi~ sex).
9.6 Kiểm định Wilcoxon cho hai mẫu (wilcox.test)
Kiểm định t dựa vào giả thiết là phân phối của một biến phải tuân theo luật phân
phối chuẩn. Nếu giả định này không đúng, kết quả của kiểm định t có thể không hợp lí
(valid). Để kiểm định phân phối của
igfi
, chúng ta có thể dùng hàm
shapiro.test
như sau:
> shapiro.test(igfi)
Shapiro-Wilk normality test
data: igfi
W = 0.8528, p-value = 1.504e-08
Trị số p nhỏ hơn 0.05 rất nhiều, cho nên chúng ta có thể nói rằng phân phối của
igfi
không tuân theo luật phân phối chuẩn. Trong trường hợp này, việc so sánh giữa hai
nhóm có thể dựa vào phương pháp phi tham số (non-parametric) có tên là kiểm định
Không có nhận xét nào:
Đăng nhận xét