Buổi 3 - Phân tích định tính

---

Giới thiệu package

Dữ liệu penguins được cung cấp bởi Tiến sĩ Kristen Gorman và Trạm Palmer, Nam Cực. Dữ liệu thu thập các đặc điểm của các loài chim cánh cụt sống ở Nam Cực. Dữ liệu này được tích hợp trong package “palmerpenguins” bao gồm 8 biến:

— species: các loại chim cánh cụt

— island : hòn đảo chim cánh cụt sinh sống

— bill_length_mm: chiều dài mỏ chim cánh cụt

— bill_depth_mm: độ sâu mỏ chim cánh cụt

— body_mass_g: khối lượng cơ thể chim cánh cụt

— sex: giống (gồm đực và cái)

— flipper_length_mm: độ dài cánh của chim cánh cụt

— year: năm thu thập số liệu

library(palmerpenguins)
d <- na.omit(penguins)
d <- as.data.frame(d)
str(d)

## 'data.frame':    333 obs. of  8 variables:
##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ bill_length_mm   : num  39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
##  $ bill_depth_mm    : num  18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
##  $ flipper_length_mm: int  181 186 195 193 190 181 195 182 191 198 ...
##  $ body_mass_g      : int  3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
##  $ year             : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
##  - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
##   ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...

Ta thấy có 3 biến định đính là Species, Island, Sex

1.Lập bảng tần số cho các biến độc lập

Biến giới tính (sex)

table(d$sex)

## 
## female   male 
##    165    168

table(d$sex)/sum(table(d$sex))

## 
##    female      male 
## 0.4954955 0.5045045

Ta thấy tỷ lệ con đực và con cái xêm nhau, gồm 168 con đực, chiếm 50.45%; 165 con cái chiếm 49.55%

Biến đảo sinh sống

table(d$island)

## 
##    Biscoe     Dream Torgersen 
##       163       123        47

table(d$island)/sum(table(d$island))

## 
##    Biscoe     Dream Torgersen 
## 0.4894895 0.3693694 0.1411411

Biscoe có tỷ lệ chim cánh cụt cao nhất (48.95%), điều này cho thấy rằng đảo Biscoe có nhiều chim cánh cụt hơn so với hai đảo còn lại. Dream có tỷ lệ chim cánh cụt trung bình (36.94%), cho thấy rằng số lượng chim cánh cụt ở đảo Dream cũng khá cao nhưng không bằng Biscoe. Torgersen có tỷ lệ chim cánh cụt thấp nhất (14.11%), cho thấy rằng đảo này có ít chim cánh cụt hơn so với hai đảo còn lại.

Biến loài

table(d$species)

## 
##    Adelie Chinstrap    Gentoo 
##       146        68       119

table(d$species)/sum(table(d$species))

## 
##    Adelie Chinstrap    Gentoo 
## 0.4384384 0.2042042 0.3573574

• Loài Adelie là phổ biến nhất trong dữ liệu, chiếm gần 44%.

• Loài Gentoo cũng khá phổ biến, chiếm gần 36%.

• Loài Chinstrap ít phổ biến nhất, chiếm khoảng 20%.

Bảng tần số giữa 2 biến

giới tính và đảo sinh sống

c <- table(d$sex,d$island)
c

##         
##          Biscoe Dream Torgersen
##   female     80    61        24
##   male       83    62        23

prop.table(c)

##         
##              Biscoe      Dream  Torgersen
##   female 0.24024024 0.18318318 0.07207207
##   male   0.24924925 0.18618619 0.06906907

— Nhận xet:

• Biscoe có tổng số chim cánh cụt nhiều nhất (163), tiếp theo là Dream (123), và ít nhất là Torgersen (47).

• Dựa vào bảng tần suất, ta thấy Số lượng chim cánh cụt đực và cái trên mỗi đảo khá cân bằng, không có sự chênh lệch lớn giữa hai giới tính.

• Biscoe có số lượng chim cánh cụt đực và cái đều cao nhất, cho thấy rằng đây có thể là môi trường sống ưa thích của chim cánh cụt trong quần đảo Palmer.

Loài và đảo sinh sống

h <- table(d$species,d$island)
h

##            
##             Biscoe Dream Torgersen
##   Adelie        44    55        47
##   Chinstrap      0    68         0
##   Gentoo       119     0         0

prop.table(h)

##            
##                Biscoe     Dream Torgersen
##   Adelie    0.1321321 0.1651652 0.1411411
##   Chinstrap 0.0000000 0.2042042 0.0000000
##   Gentoo    0.3573574 0.0000000 0.0000000

• Nhận xét:

• Trên đảo Biscoe, loài Gentoo chiếm tỷ lệ cao nhất với 35.74%, trong khi loài Adelie chiếm 13.21%, và không có loài Chinstrap nào.

• Trên đảo Dream, loài Chinstrap chiếm tỷ lệ cao nhất với 20.42%, theo sau là loài Adelie với 16.52%, và không có loài Gentoo nào.

• Trên đảo Torgersen, chỉ có loài Adelie xuất hiện với tỷ lệ 14.11%, và không có loài Chinstrap hay Gentoo nào.

2. Vẽ đồ thị

2.1 Đồ thị cột

Vẽ đồ thị cột cho Biến loài (species)

library(ggplot2)
ggplot(d, aes(species,fill=species)) + geom_bar(position = 'dodge')

Nhìn vào biểu đồ ta thấy loài Adelie chiếm tỷ lệ nhiều nhất, kế đến là Gentoo và cuối cùng là Chinstrap.

Chúng ta có sắp xếp lại biểu đồ theo thứ tự tăng dần và hiển thị số cá thể trên từng loài:

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

d %>%
  filter(species %in% c("Adelie", "Chinstrap", "Gentoo")) %>%
  group_by(species) %>%
  count() %>%
  ggplot(aes(reorder(species,-n), n)) +
  geom_col(width = 0.7,fill="brown") +
  geom_text(aes(label = n), vjust = 2, color = "white", size = 4) +
  labs(x = NULL, y = "Number of Penguins",
       title = "Number of Penguins by Species",
       caption = "Data Source: Custom Penguins Dataset") ->> p1

# Hiển thị biểu đồ
p1

Quay ngược biểu đồ, như sau:

p1 + coord_flip()

Lúc này với biểu đồ cột ngang, các số có trên biểu đồ bị mất. Do đó ta hiệu chỉnh như sau:

Hiệu chỉnh như sau:

penguins_count <- d %>%
  group_by(species) %>%
  count()
theme_set(theme_minimal())
p2 <- ggplot(penguins_count, aes(reorder(species, n), n)) +
  geom_col(width = 0.7,fill="brown") +
  geom_text(aes(label = n), hjust = 1.2, color = "white", size = 5) +
  labs(x = NULL, y = "Number of Penguins",
       title = "Number of Penguins by Species",
       caption = "Data Source: Custom Penguins Dataset") +
  coord_flip() +
  scale_y_continuous(breaks = seq(0, max(penguins_count$n), 10))
p2

Sử dụng theme của một số tạp chí nổi tiếng:

library(ggthemes)
p2 + 
  theme_economist(horizontal = FALSE)+
  scale_fill_economist()

library(dplyr)

d %>%
  filter(species %in% c("Adelie", "Chinstrap", "Gentoo")) %>%
  group_by(species) %>%
  count() %>%
  ggplot(aes(reorder(species, n), n, fill = species)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  geom_text(aes(label = n), hjust = 1.2, color = "white", size = 5) +
  labs(x = NULL,
       y = "Số lượng",
       title = "Số lượng chim cánh cụt theo loài",
       caption = "Nguồn dữ liệu: https://github.com/allisonhorst/palmerpenguins") +
  coord_flip() +
  scale_y_continuous(breaks = seq(0, 100, 10)) -> p3
print(p3)

library(ggthemes)
p3 + 
  theme_economist(horizontal = FALSE) + 
  scale_fill_economist()

— Nhận xét: Đồ thị đã được cải tiến một cách hoàn chỉnh và đẹp hơn. Nhìn vào biểu đồ ta thấy loài Adelie chiếm tỷ lệ nhiều nhất, với 146 cá thể kế đến là Gentoo với 119 cá thể và cuối cùng là Chinstrap với 68 cá thể.

3. Ước lượng tỷ lệ

— Công thức ước lượng tỉ lệ (cho 1 tổng thể):

\[ \hat{p} - Z_{\alpha/2} \sqrt{\frac{\hat{p}(1 - \hat{p})}{n}} \le P \le \hat{p} + Z_{\alpha/2} \sqrt{\frac{\hat{p}(1 - \hat{p})}{n}} \]

Trong đó:

• \(\hat{p}\) là tỷ lệ mẫu của mức độ ngoại hình.
• \(P\) là tỷ lệ số lượng của từng phân loại ngoại hình.
• \(Z_{\alpha/2}\) là giá trị critical từ phân phối chuẩn tương ứng với mức tin cậy \(1 - \alpha\).
• \(n\) là kích thước mẫu.

Ước lượng giống cánh cụt đực, đồng thời kiểm định xem tỷ lệ % cá thể có đạt 50% hay không (giả thuyết Ho=50%)

mal <- d[d$sex == "male",]
prop.test(length(d$sex), length(d$sex), p = 0.5)

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(d$sex) out of length(d$sex), null probability 0.5
## X-squared = 331, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.9857837 1.0000000
## sample estimates:
## p 
## 1

Vì p_value = 2.2e-16 < 1%. Do đó, bác bỏ H0. Vậy tỷ lệ con đực không đạt 50%.

Ước lượng các loài sinh sống ở Đảo Biscoe, đồng thời kiểm định xem số cá thể sống ở Biscoe có đạt 40% hay không (Giả thuyết H0=0.4

Ade <- d[d$island == "Biscoe",]
prop.test(length(d$island), length(d$island), p = 0.4)

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(d$island) out of length(d$island), null probability 0.4
## X-squared = 497, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.4
## 95 percent confidence interval:
##  0.9857837 1.0000000
## sample estimates:
## p 
## 1

Vì p_value = 2.2e-16 < 1%. Do đó, bác bỏ H0. Vậy tỷ lệ cá thể sống ở đảo Bisoce không đạt 40%.

4.Ước lượng chênh lệch 2 tỷ lệ

Trong đó:

• \(\hat{p}_1\) và \(\hat{p}_2\) là tỷ lệ chỉ số xếp hạng ngoại hình cao hơn 1 trong mẫu nam và nữ tương ứng.
• \(p_1\) và \(p_2\) là tỷ lệ chỉ số xếp hạng ngoại hình cao hơn 1 trong toàn bộ tổng thể nam và nữ tương ứng.
• \(n_1\) và \(n_2\) là kích thước của mẫu nam và nữ tương ứng.
• \(Z_{\alpha/2}\) là giá trị critical từ phân phối chuẩn tương ứng với mức tin cậy \(1 - \alpha\).

Thực hiện bài toán kiểm định giả thuyết sự bằng nhau về tỷ lệ cá thể đực và cái sống ở đảo Bisoce 2 tổng thể, nghĩa là chúng ta thực hiện bài toán kiểm định H0:p1=p2

d4 <- d[d$sex == "male",]
d4f <- d[d$sex=="female",]
d4d4 <- d4[d4$island =="Biscoe",]
d4fd4 <- d4f[d4f$island == "Biscoe",]

a <- c(nrow(d4), nrow(d4f))
b <- c(nrow(d4d4), nrow(d4fd4))

prop.test(b,a)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  b out of a
## X-squared = 0.0033955, df = 1, p-value = 0.9535
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1041884  0.1225867
## sample estimates:
##    prop 1    prop 2 
## 0.4940476 0.4848485

với p_value = 0.9535, chấp nhận H0. Vậy tỷ lệ cá thể đực ở đảo Biscoe bằng với tỷ lệ cái sống ở đảo Biscoe.

5. Tính rủi ro tương đối

— Trong phần này, tôi tiến hành tính rủi ro tương đối

h <- table(d$species,d$island)
addmargins(h)

##            
##             Biscoe Dream Torgersen Sum
##   Adelie        44    55        47 146
##   Chinstrap      0    68         0  68
##   Gentoo       119     0         0 119
##   Sum          163   123        47 333

Vì bảng tần số có dạng ma trận 3x3, do đó ta không thể dùng Relrisk để tính toán rủi ro tương đối được, do đó ta tính thủ công từng cặp với nhau như sau:

Cố định loài Adelie, ta sẽ tính rủi ro tương đối của loài Adelie so với 3 đảo. Dưới đây là dữ liệu từ bảng tần số:

# Dữ liệu
data <- data.frame(
  Species = c("Adelie", "Chinstrap", "Gentoo"),
  Biscoe = c(44, 0, 119),
  Dream = c(55, 68, 0),
  Torgersen = c(47, 0, 0),
  Sum = c(146, 68, 119)
)
data

Thực hiện tính toán, như sau:

P_Biscoe_Adelie <- data$Biscoe[1] / data$Sum[1]
P_Dream_Adelie <- data$Dream[1] / data$Sum[1]
P_Torgersen_Adelie <- data$Torgersen[1] / data$Sum[1]
P_Biscoe_Adelie

## [1] 0.3013699

P_Dream_Adelie

## [1] 0.3767123

P_Torgersen_Adelie

## [1] 0.3219178

Rủi ro tương đối cho Adelie

# Rủi ro tương đối cho Adelie
RR_1 <- P_Biscoe_Adelie / P_Dream_Adelie
RR_1

## [1] 0.8

• Kết luận: Rủi ro tương đối cho Adelie (RR_1) là khoảng 0.80. Điều này có nghĩa là Adelie ở đảo Biscoe có nguy cơ phân bố chỉ bằng 0.80 lần so với Adelie ở đảo Dream.

— Tương tự, ta có thể tính rủi ro tương đối của loài Adelie ở Đảo Biscoe so với Đảo Torgersen như sau:

RR_2 <- P_Biscoe_Adelie / P_Torgersen_Adelie
RR_2

## [1] 0.9361702

• Kết luận: Rủi ro tương đối cho Adelie (RR_2) là khoảng 0.94. Điều này có nghĩa là Adelie ở đảo Biscoe có nguy cơ phân bố chỉ bằng 0.94 lần so với Adelie ở đảo Torgersen.

— Tương tự, ta có thể tính rủi ro tương đối của loài Adelie ở Đảo Dream so với Đảo Torgersen như sau:

RR_3 <- P_Dream_Adelie / P_Torgersen_Adelie
RR_3

## [1] 1.170213

• Kết luận: Rủi ro tương đối cho Adelie (RR_3) là khoảng 1.17. Điều này có nghĩa là Adelie ở đảo Dream có nguy cơ phân bố bằng 1.17 lần so với Adelie ở đảo Dream.

6. Ước lượng odd ratio

Đầu tiên ta lập bảng tần số cho 2 biến sex và species

odd <- table(d$species,d$sex)
addmargins(odd)

##            
##             female male Sum
##   Adelie        73   73 146
##   Chinstrap     34   34  68
##   Gentoo        58   61 119
##   Sum          165  168 333

Tương tự, vì đây là ma trận 3x2 do đó ta không dùng câu lệnh để tính odd ratio được. Thay vào đó, ta làm thủ công cho từng cặp biểu hiện của từng biến với nhau.

Làm mới lại bảng tần số để tính.

da6 <- matrix(c(73, 73, 34, 34, 58, 61), nrow = 3, byrow = TRUE,
               dimnames = list(c("Adelie", "Chinstrap", "Gentoo"),
                               c("female", "male")))
rownames(data) <- c("Adelie", "Chinstrap", "Gentoo")
colnames(data) <- c("female", "male")
da6

##           female male
## Adelie        73   73
## Chinstrap     34   34
## Gentoo        58   61

Tính Odds-ratio tỷ lệ nam nữ của loài Adelie và Chintrap xem liệu rằng tỷ lệ nam/ nữ của 2 loài này có sự khác biệt hay không

a_adelie <- da6["Adelie", "female"]
b_adelie <- da6["Adelie", "male"]
c_chinstrap <- da6["Chinstrap", "female"]
d_chinstrap <- da6["Chinstrap", "male"]

OR_ac <- (a_adelie / b_adelie) / (c_chinstrap / d_chinstrap)
OR_ac

## [1] 1

— Nhận xét: Vậy với odds ratio =1 thì tỷ lệ nam nữ ở cả 2 loài Adelie và Chinstrap là như nhau.

Tương tự, ta có thể tính odds ratio cho loài Chinstrap và Gentoo như sau:

e_gentoo <- da6["Gentoo","female"]
f_gentoo <- da6["Gentoo","male"]
OR_cg <- (c_chinstrap / d_chinstrap)/(e_gentoo/f_gentoo)
OR_cg

## [1] 1.051724

— Nhận xét: Khi tính toán Odds Ratio giữa loài Chinstrap và Gentoo và kết quả là 1.051724, điều này có nghĩa tồn tại sự chênh lệch giữa cá thể đực và cái của cả 2 loài. Cụ thể, sự chênh lệch đến từ việc loài Gentoo có 61 cá thể đực nhiều hơn cá thể cái (58 cá thể).