#1.2. Mô tả đặc điểm tuổi (age), giới tính  (gender), cân nặng (weight), chiều cao (height), tỉ trọng mỡ  (pcfat), tiền căn bệnh cao huyết áp (hypertension)và tiền căn bệnh tiểu đường  (diabetes)
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
table1(~ age + gender + weight + height + pcfat + hypertension + diabetes, data = ob)
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Min, Max] 48.0 [13.0, 88.0]
gender
F 862 (70.8%)
M 355 (29.2%)
weight
Mean (SD) 55.1 (9.40)
Median [Min, Max] 54.0 [34.0, 95.0]
height
Mean (SD) 157 (7.98)
Median [Min, Max] 155 [136, 185]
pcfat
Mean (SD) 31.6 (7.18)
Median [Min, Max] 32.4 [9.20, 48.4]
hypertension
Mean (SD) 0.507 (0.500)
Median [Min, Max] 1.00 [0, 1.00]
diabetes
Mean (SD) 0.111 (0.314)
Median [Min, Max] 0 [0, 1.00] 
#1.3 Bạn nhận xét như thế nào về kết quả của tiền căn bệnh cao huyết áp và tiểu đường. Làm cách nào để trình bày kết quả tốt hơn?
ob$hyper = as.factor(ob$hypertension)
ob$diab = as.factor(ob$diabetes)

table1(~ age + gender + weight + height + pcfat + hypertension + hyper + diabetes + diab, data = ob)
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Min, Max] 48.0 [13.0, 88.0]
gender
F 862 (70.8%)
M 355 (29.2%)
weight
Mean (SD) 55.1 (9.40)
Median [Min, Max] 54.0 [34.0, 95.0]
height
Mean (SD) 157 (7.98)
Median [Min, Max] 155 [136, 185]
pcfat
Mean (SD) 31.6 (7.18)
Median [Min, Max] 32.4 [9.20, 48.4]
hypertension
Mean (SD) 0.507 (0.500)
Median [Min, Max] 1.00 [0, 1.00]
hyper
0 600 (49.3%)
1 617 (50.7%)
diabetes
Mean (SD) 0.111 (0.314)
Median [Min, Max] 0 [0, 1.00]
diab
0 1082 (88.9%)
1 135 (11.1%) 
# Tỷ lệ người có tiền căn bệnh cao huyết áp chiếm 50.7%, điều này cho thấy đây là một bệnh phổ biến trong mẫu nghiên cứu này
# Khoảng 11% người tham gia có tiền căn tiểu đường, tức là tỷ lệ cao hơn một chút so với người bị cao huyết áp trong cùng mẫu dữ liệu phân tích
# Cách trình bày tốt hơn: thay vì đưa dữ liệu thô có thể trình bày qua bảng tóm tắt hoặc qua biểu đồ
library(ggplot2)

# Gender
ggplot(ob, aes(x = gender)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat='count', aes(label=..count..), vjust=-0.5) +
  labs(title="Phân bố giới tính", x="Giới tính", y="Số lượng") +
  theme_minimal()
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Hypertension
ggplot(ob, aes(x = factor(hypertension, labels=c("No","Yes")))) +
  geom_bar(fill = "tomato") +
  geom_text(stat='count', aes(label=..count..), vjust=-0.5) +
  labs(title="Tăng huyết áp", x="", y="Số lượng") +
  theme_minimal()

# Boxplot age theo giới tính
ggplot(ob, aes(x = gender, y = age, fill = gender)) +
  geom_boxplot() +
  labs(title="Tuổi theo giới tính", x="Giới tính", y="Tuổi") +
  theme_minimal()

# Boxplot pcfat (tỷ lệ mỡ cơ thể) theo giới tính
ggplot(ob, aes(x = gender, y = pcfat, fill = gender)) +
  geom_boxplot() +
  labs(title="Tỷ lệ mỡ cơ thể theo giới tính", x="Giới tính", y="% mỡ cơ thể") +
  theme_minimal()

#1.4.  Trình bày kết quả trung vị (Q1, Q3) thay vi trung vị (min, max) cho biến liên tục.
table1(~ age + weight + height + pcfat, data = ob, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]"))
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Q1, Q3] 48.0 [35.0, 58.0]
weight
Mean (SD) 55.1 (9.40)
Median [Q1, Q3] 54.0 [49.0, 61.0]
height
Mean (SD) 157 (7.98)
Median [Q1, Q3] 155 [151, 162]
pcfat
Mean (SD) 31.6 (7.18)
Median [Q1, Q3] 32.4 [27.0, 36.8] 
#1.5. Mô tả theo giới tính
table1(~ age + weight + height + pcfat + hyper + diab | gender, data = ob)
F
(N=862)		 M
(N=355)		 Overall
(N=1217)
age
Mean (SD) 48.6 (16.4) 43.7 (18.8) 47.2 (17.3)
Median [Min, Max] 49.0 [14.0, 85.0] 44.0 [13.0, 88.0] 48.0 [13.0, 88.0]
weight
Mean (SD) 52.3 (7.72) 62.0 (9.59) 55.1 (9.40)
Median [Min, Max] 51.0 [34.0, 95.0] 62.0 [38.0, 95.0] 54.0 [34.0, 95.0]
height
Mean (SD) 153 (5.55) 165 (6.73) 157 (7.98)
Median [Min, Max] 153 [136, 170] 165 [146, 185] 155 [136, 185]
pcfat
Mean (SD) 34.7 (5.19) 24.2 (5.76) 31.6 (7.18)
Median [Min, Max] 34.7 [14.6, 48.4] 24.6 [9.20, 39.0] 32.4 [9.20, 48.4]
hyper
0 430 (49.9%) 170 (47.9%) 600 (49.3%)
1 432 (50.1%) 185 (52.1%) 617 (50.7%)
diab
0 760 (88.2%) 322 (90.7%) 1082 (88.9%)
1 102 (11.8%) 33 (9.3%) 135 (11.1%) 
#1.6. Đánh giá xem đặc điểm nào là khác biệt đáng kể (significant difference) giữa nam và nữ.
library(compareGroups)
createTable(compareGroups(gender ~ age + weight + height + pcfat + hyper + diab, data = ob))
## 
## --------Summary descriptives table by 'gender'---------
## 
## ________________________________________ 
##             F           M      p.overall 
##           N=862       N=355              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## age    48.6 (16.4) 43.7 (18.8)  <0.001   
## weight 52.3 (7.72) 62.0 (9.59)  <0.001   
## height 153 (5.55)  165 (6.73)   <0.001   
## pcfat  34.7 (5.19) 24.2 (5.76)  <0.001   
## hyper:                           0.569   
##     0  430 (49.9%) 170 (47.9%)           
##     1  432 (50.1%) 185 (52.1%)           
## diab:                            0.238   
##     0  760 (88.2%) 322 (90.7%)           
##     1  102 (11.8%) 33 (9.30%)            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

việc 2: Phân tích khác biệt giữa 2 nhóm

#2,1 Dữ liệu về tải trọng của 2 nhóm A và B như sau
A = c(14, 4, 10, 6, 3, 11, 12)
B = c(16, 17, 13, 12, 7, 16, 11, 8, 7)
wt = c(A, B)
group = c(rep("A", 7), rep("B", 9))
df = data.frame(wt, group)
dim(df)
## [1] 16  2
#2.2. Tải trọng có tuân theo phân bổ chuẩn (normal distribution) không?
library(lessR)
## 
## lessR 4.4.5                         feedback: gerbing@pdx.edu 
## --------------------------------------------------------------
## > d <- Read("")  Read data file, many formats available, e.g., Excel
##   d is default data frame, data= in analysis routines optional
## 
## Many examples of reading, writing, and manipulating data, 
## graphics, testing means and proportions, regression, factor analysis,
## customization, forecasting, and aggregation from pivot tables
##   Enter: browseVignettes("lessR")
## 
## View lessR updates, now including time series forecasting
##   Enter: news(package="lessR")
## 
## Interactive data analysis
##   Enter: interact()
## 
## Attaching package: 'lessR'
## The following object is masked from 'package:table1':
## 
##     label
Histogram(wt, data = df)
## >>> Note: wt is not in a data frame (table)
## >>> Note: wt is not in a data frame (table)

## >>> Suggestions 
## bin_width: set the width of each bin 
## bin_start: set the start of the first bin 
## bin_end: set the end of the last bin 
## Histogram(wt, density=TRUE)  # smoothed curve + histogram 
## Plot(wt)  # Violin/Box/Scatterplot (VBS) plot 
## 
## --- wt --- 
##  
##      n   miss     mean       sd      min      mdn      max 
##      16      0    10.44     4.29     3.00    11.00    17.00 
##  
## 
## No (Box plot) outliers 
## 
## 
## Bin Width: 2 
## Number of Bins: 8 
##  
##      Bin  Midpnt  Count    Prop  Cumul.c  Cumul.p 
## ------------------------------------------------- 
##   2 >  4       3      2    0.12        2     0.12 
##   4 >  6       5      1    0.06        3     0.19 
##   6 >  8       7      3    0.19        6     0.38 
##   8 > 10       9      1    0.06        7     0.44 
##  10 > 12      11      4    0.25       11     0.69 
##  12 > 14      13      2    0.12       13     0.81 
##  14 > 16      15      2    0.12       15     0.94 
##  16 > 18      17      1    0.06       16     1.00 
## 
#Kiểm định test
shapiro.test(df$wt)
## 
##  Shapiro-Wilk normality test
## 
## data:  df$wt
## W = 0.96213, p-value = 0.7006
# Kết quả kiểm định test cho kết quả W = 0.96213 chỉ số Shapiro - Wilk càng gần 1 thì dữ liệu càng giống phân phối chuẩn,  p -value = 0.7006>0.05 tức xác suất quan sát được thống kê W nếu dữ liệu thực sự phân phối chuẩn, không bác bỏ giả thiết H0
# Histogram
hist(df$wt, main="Histogram of wt", xlab="wt", col="lightblue", border="black")

# Q-Q plot
qqnorm(df$wt)       # vẽ Q-Q plot
qqline(df$wt, col="red")  # thêm đường chuẩn

#2.3. Mô tả đặc điểm tải trọng
library(table1)
table1(~ wt | group, data = df, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]"))
A
(N=7)		 B
(N=9)		 Overall
(N=16)
wt
Mean (SD) 8.57 (4.24) 11.9 (3.95) 10.4 (4.29)
Median [Q1, Q3] 10.0 [5.00, 11.5] 12.0 [8.00, 16.0] 11.0 [7.00, 13.3] 
#Kết quả cho thấy nhóm B có xu hướng wt cao hơn nhóm A. Cần kiểm định thống kê t để xác định sự khác biệt có ý nghĩa hay không
#2.4. Thực hiện phép kiểm định t
t.test(A, B)
## 
##  Welch Two Sample t-test
## 
## data:  A and B
## t = -1.6, df = 12.554, p-value = 0.1345
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.813114  1.178194
## sample estimates:
## mean of x mean of y 
##  8.571429 11.888889
t.test(wt ~ group, data = df)
## 
##  Welch Two Sample t-test
## 
## data:  wt by group
## t = -1.6, df = 12.554, p-value = 0.1345
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
##  -7.813114  1.178194
## sample estimates:
## mean in group A mean in group B 
##        8.571429       11.888889
# Kết quả p-value = 0.1345 > 0.05, không bác bỏ H0. Tức không có sự khác biệt có ý nghĩa thống kê giữa trung bình wt của 2 nhóm ở mức ý nghĩa 5%. Đồng thời khoảng tin cậy [-7.81, 1.18] cũng bao gồm 0 như vậy củng cố thêm kết luận trên
#2.5. Thực hiện bootstrap
library(simpleboot)
## Simple Bootstrap Routines (1.1-8)
library(boot)
b = two.boot(A, B, mean, R = 1000)
boot.ci(b)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-7.222,  0.367 )   (-7.633,  0.317 )  
## 
## Level     Percentile            BCa          
## 95%   (-6.952,  0.998 )   (-6.947,  1.015 )  
## Calculations and Intervals on Original Scale
hist(b, breaks = 50)

#2.6. Sử dụng chatGPT
#Câu lệnh prombt: Tôi có dữ liệu cho 2 nhóm như sau. Bạn có thể viết mã R để phân tích sự khác biệt giữa hai nhóm dùng bootstrap? A = c(14, 4, 10, 6, 3, 11, 12) B = c(16, 17, 13, 12, 7, 16, 11, 8, 7
library(boot)

# Dữ liệu
A <- c(14, 4, 10, 6, 3, 11, 12)
B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)

# Kết hợp thành data frame
df <- data.frame(
  value = c(A, B),
  group = factor(c(rep("A", length(A)), rep("B", length(B))))
)

# Hàm tính khác biệt trung bình giữa 2 nhóm
diff_mean <- function(data, indices) {
  d <- data[indices, ]  # lấy mẫu bootstrap
  mean_A <- mean(d$value[d$group == "A"])
  mean_B <- mean(d$value[d$group == "B"])
  return(mean_A - mean_B)
}

# Thực hiện bootstrap với 5000 lần lặp
set.seed(123)  # để kết quả có thể lặp lại
b <- boot(data = df, statistic = diff_mean, R = 5000)

# Xem kết quả bootstrap
print(b)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = df, statistic = diff_mean, R = 5000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* -3.31746 -0.02336054    1.991425
# Khoảng tin cậy 95% theo phương pháp percentile
boot.ci(b, type = "perc")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 4999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   (-7.205,  0.567 )  
## Calculations and Intervals on Original Scale