18  Các kiểm định thống kê cơ bản

Chương này sẽ trình bày cách để thực hiện các phép kiểm định thống kê cơ bản bằng cách sử dụng base R, rstatix, và gtsummary.

…nhiều kiểm định khác có thể được thực hiện, nhưng chúng tôi chỉ trình bày các kiểm định thông dụng và kết nối với các phần khác trong cuốn sổ tay này.

Mỗi package được đề cập bên trên đều có một số ưu điểm và khuyết điểm nhất định:

18.1 Các bước chuẩn bị

Gọi các packages

Đoạn code này hiển thị việc gọi các package cần thiết cho phân tích. Trong cuốn sổ tay này, chúng tôi nhấn mạnh đến hàm p_load() trong package pacman, cài đặt gói lệnh nếu cần thiết gọi chúng ra để sử dụng. Các package đã cài đặt cũng có thể được gọi ra bằng library() từ base R. Xem thêm thông tin các package của R trong chương R cơ bản.

pacman::p_load(
  rio,          # File import
  here,         # File locator
  skimr,        # get overview of data
  tidyverse,    # data management + ggplot2 graphics, 
  gtsummary,    # summary statistics and tests
  rstatix,      # statistics
  corrr,        # correlation analayis for numeric variables
  janitor,      # adding totals and percents to tables
  flextable     # converting tables to HTML
  )

Nhập số liệu

Chúng ta nhập bộ số liệu của các ca bệnh về một vụ dịch Ebola mô phỏng. Để tiện theo dõi, bấm để tải bộ số liệu linelist “đã được làm sạch” (as .rds file). Nhập số liệu bằng hàm import() từ package rio package (nó chấp nhận nhiều loại tập tin như .xlsx, .rds, .csv - xem thêm chương Nhập xuất dữ liệu để biết thêm chi tiết).

# import the linelist
linelist <- import("linelist_cleaned.rds")

50 hàng đầu tiên của bộ dữ liệu linelist được hiển thị như dưới đây.

18.2 Các kiểm định trong base R

Các lệnh trong base R functions to conduct statistical tests. có thể được sử dụng để thực hiện các kiểm định thống kê. Các câu lệnh tương đối đơn giản và kết quả sẽ hiển thị trong bảng điều khiển R Console. Tuy nhiên, kết quả đầu ra thường dưới dạng liệt kê, vì thế sẽ khó thao tác hơn nếu muốn sử dụng kết quả trong các thao tác tiếp theo.

Kiểm định t

Một kiểm định t, hay còn được gọi là “Student’s t-Test”, thường được sử dụng để xác định có sự khác biệt có ý nghĩa thống kê giữa giá trị trung bình của hai nhóm. Bên dưới là cú pháp để thực hiện kiểm định này tùy thuộc vào các cột có trong cùng một data frame hay không.

Cú pháp 1: Đây là cú pháp khi cột của biến liên tục và phân loại nằm trong cùng một data frame. Đặt biến liên tục bên trái và biến phân loại bên phải của phương trình. Ghi rõ bộ số liệu sau data =. Các tùy chọn khác như số liệu bắt cặp, viết thêm paired = TRUE, khoảng tin cậy, viết thêm conf.level = (mặc định là 0.95), và giả thuyết thay thế alternative = (hai đuôi - “two.sided”, hoặc một đuôi nhỏ hơn hay lớn hơn - “less”, or “greater”). Gõ ?t.test để biết thêm chi tiết.

## compare mean age by outcome group with a t-test
t.test(age_years ~ gender, data = linelist)

    Welch Two Sample t-test

data:  age_years by gender
t = -21.344, df = 4902.3, p-value < 2.2e-16
alternative hypothesis: true difference in means between group f and group m is not equal to 0
95 percent confidence interval:
 -7.571920 -6.297975
sample estimates:
mean in group f mean in group m 
       12.60207        19.53701 

Cú pháp 2: Đây là cú pháp khi so sánh hai véc tơ dạng số. Ví dụ như hai cột nằm trong hai bộ số liệu khác nhau.

t.test(df1$age_years, df2$age_years)

Kiểm định t cũng được sử dụng để xác định có sự khác biệt có ý nghĩa thống kê giữa giá trị trung bình của mẫu với một số giá trị cụ thể. Đây là phép kiểm định t cho một mẫu với trung bình quần thể giả thuyết/đã biết như mu =:

t.test(linelist$age_years, mu = 45)

Kiểm định Shapiro-Wilk

Kiểm định Shapiro-Wilk có thể được sử để xác định xem một mẫu có phân bố bình thường/phân bố chuản hay không (một giả định của nhiều kiểm định khác, ví dụ như kiểm định t). Tuy nhiên, phép kiểm định này chỉ có thể được sử dụng cho một mẫu có từ 3 đến 5000 quan sát. Đối với cỡ mẫu lớn hơn, nên sử dụng biểu đồ quantile-quantile plot.

shapiro.test(linelist$age_years)

Kiểm định tổng thứ hạng Wilcoxon

Kiểm định tổng thứ hạng Wilcoxon, hay còn gọi là kiểm định Mann–Whitney U, thường được sử dụng để giúp xác định xem hai mẫu có cùng phân bố hay không khi quần thể của chúng không có phân bố chuẩn hoặc có phương sai không bằng nhau.

## compare age distribution by outcome group with a wilcox test
wilcox.test(age_years ~ outcome, data = linelist)

    Wilcoxon rank sum test with continuity correction

data:  age_years by outcome
W = 2501868, p-value = 0.8308
alternative hypothesis: true location shift is not equal to 0

Kiểm định Kruskal-Wallis

Kiểm định Kruskal-Wallis là một phần mở rộng của kiểm định tổng thứ hạng Wilcoxon mà có thể được sử dụng để kiểm định sự khác biệt trong phân bố của nhiều hơn hai mẫu. Khi có hai mẫu được sử dụng, nó cho kết quả giống như của kiểm định tổng thứ hạng Wilcoxon.

## compare age distribution by outcome group with a kruskal-wallis test
kruskal.test(age_years ~ outcome, linelist)

    Kruskal-Wallis rank sum test

data:  age_years by outcome
Kruskal-Wallis chi-squared = 0.045675, df = 1, p-value = 0.8308

Kiểm định Chi bình phương

Kiểm định Chi bình phương của Pearson được sử dụng trong kiểm tra sự khác biệt có ý nghĩa thống kê giữa các biến phân loại.

## compare the proportions in each group with a chi-squared test
chisq.test(linelist$gender, linelist$outcome)

    Pearson's Chi-squared test with Yates' continuity correction

data:  linelist$gender and linelist$outcome
X-squared = 0.0011841, df = 1, p-value = 0.9725

18.3 rstatix package

Package rstatix cho phép thực hiện các kiểm định thống kê và truy xuất kết quả “dễ sử dụng cho các tính toán tiếp theo”. Có nghĩa là kết quả xuất tự động thành một data frame để có thể thực hiện các thao tác tiếp theo. Nó cũng dễ dàng để nhóm dữ liệu mà sẽ được chuyền vào các hàm, ở đó các thống kê được thực hiện cho từng nhóm.

Tóm tắt thống kê

Hàm get_summary_stats() là một cách thực hiện tóm tắt thống kê nhanh. Chỉ cần đưa bộ số liệu và chỉ định các cột muốn phân tích vào hàm này. Nếu không có cột nào được cụ thể, tóm tắt thống kê sẽ tính toán cho tất cả các cột.

Tóm tắt thống kê đầy đủ sẽ cho kết quả mặc định như sau: số quan sát (n), giá trị nhỏ nhất, giá trị lớn nhất, trung vị, giá trị tứ phân vị thứ nhất (25%), giá trị tứ phân vị thứ ba (75%), khoảng tứ phân vị, độ lệch tuyệt đối của trung vị (mad), trung bình, độ lệch chuẩn, sai số chuẩn và khoảng tin cậy của trung bình.

linelist %>%
  rstatix::get_summary_stats(age, temp)
# A tibble: 2 × 13
  variable     n   min   max median    q1    q3   iqr    mad  mean     sd    se
  <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>
1 age       5802   0    84     13     6    23      17 11.9    16.1 12.6   0.166
2 temp      5739  35.2  40.8   38.8  38.2  39.2     1  0.741  38.6  0.977 0.013
# ℹ 1 more variable: ci <dbl>

Có thể tóm tắt một số giá trị thống kê bằng cách cung cấp một trong số các giá trị sau đến type =: “full”, “common”, “robust”, “five_number”, “mean_sd”, “mean_se”, “mean_ci”, “median_iqr”, “median_mad”, “quantile”, “mean”, “median”, “min”, “max”.

Nó cũng có thể được sử dụng để nhóm số liệu, sao cho một hàng được trả về cho mỗi biến nhóm:

linelist %>%
  group_by(hospital) %>%
  rstatix::get_summary_stats(age, temp, type = "common")
# A tibble: 12 × 11
   hospital     variable     n   min   max median   iqr  mean     sd    se    ci
   <chr>        <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
 1 Central Hos… age        445   0    58     12    15    15.7 12.5   0.591 1.16 
 2 Central Hos… temp       450  35.2  40.4   38.8   1    38.5  0.964 0.045 0.089
 3 Military Ho… age        884   0    72     14    18    16.1 12.4   0.417 0.818
 4 Military Ho… temp       873  35.3  40.5   38.8   1    38.6  0.952 0.032 0.063
 5 Missing      age       1441   0    76     13    17    16.0 12.9   0.339 0.665
 6 Missing      temp      1431  35.8  40.6   38.9   1    38.6  0.97  0.026 0.05 
 7 Other        age        873   0    69     13    17    16.0 12.5   0.422 0.828
 8 Other        temp       862  35.7  40.8   38.8   1.1  38.5  1.01  0.034 0.067
 9 Port Hospit… age       1739   0    68     14    18    16.3 12.7   0.305 0.598
10 Port Hospit… temp      1713  35.5  40.6   38.8   1.1  38.6  0.981 0.024 0.046
11 St. Mark's … age        420   0    84     12    15    15.7 12.4   0.606 1.19 
12 St. Mark's … temp       410  35.9  40.6   38.8   1.1  38.5  0.983 0.049 0.095

Bạn cũng có thể sử dụng rstatix để thực hiện các kiểm định thống kê:

Kiểm định t

USử dụng cú pháp để chỉ định cột biến liên tục và cột biến phân loại:

linelist %>% 
  t_test(age_years ~ gender)
# A tibble: 1 × 10
  .y.   group1 group2    n1    n2 statistic    df        p    p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl> <chr>       
1 age_… f      m       2807  2803     -21.3 4902. 9.89e-97 9.89e-97 ****        

Hoặc sử dụng ~ 1 và ghi rõ mu = cho kiểm định t một mẫu. Cú pháp này có thể sử dụng để thực hiện cho nhóm.

linelist %>% 
  t_test(age_years ~ 1, mu = 30)
# A tibble: 1 × 7
  .y.       group1 group2         n statistic    df     p
* <chr>     <chr>  <chr>      <int>     <dbl> <dbl> <dbl>
1 age_years 1      null model  5802     -84.2  5801     0

Nếu có thể, các kiểm định thống kê có thể thực hiện theo nhóm, như được trình bày bên dưới.

linelist %>% 
  group_by(gender) %>% 
  t_test(age_years ~ 1, mu = 18)
# A tibble: 3 × 8
  gender .y.       group1 group2         n statistic    df         p
* <chr>  <chr>     <chr>  <chr>      <int>     <dbl> <dbl>     <dbl>
1 f      age_years 1      null model  2807    -29.8   2806 7.52e-170
2 m      age_years 1      null model  2803      5.70  2802 1.34e-  8
3 <NA>   age_years 1      null model   192     -3.80   191 1.96e-  4

Kiểm định Shapiro-Wilk

Như đã đề cập bên trên, cỡ mẫu phải nằm trong khoảng từ 3 đến 5000.

linelist %>% 
  head(500) %>%            # first 500 rows of case linelist, for example only
  shapiro_test(age_years)
# A tibble: 1 × 3
  variable  statistic        p
  <chr>         <dbl>    <dbl>
1 age_years     0.917 6.67e-16

Kiểm định tổng thứ hạng Wilcoxon

linelist %>% 
  wilcox_test(age_years ~ gender)
# A tibble: 1 × 9
  .y.       group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
* <chr>     <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
1 age_years f      m       2807  2803   2829274 3.47e-74 3.47e-74 ****        

Kiểm định Kruskal-Wallis

Cũng được biết như kiểm định Mann-Whitney U.

linelist %>% 
  kruskal_test(age_years ~ outcome)
# A tibble: 1 × 6
  .y.           n statistic    df     p method        
* <chr>     <int>     <dbl> <int> <dbl> <chr>         
1 age_years  5888    0.0457     1 0.831 Kruskal-Wallis

Kiểm định Chi bình phương

Hàm kiểm định Chi bình phương chấp nhận một bảng, vì vậy đầu tiên là tạo một bảng chéo. Có nhiều cách để tạo một bảng chéo (xem chương Bảng mô tả) nhưng ở đây chúng ta sử dụng hàm tabyl() từ janitor avà bỏ cột ngoài cùng bên trái của nhãn giá trị trước khi đưa vào hàm chisq_test().

linelist %>% 
  tabyl(gender, outcome) %>% 
  select(-1) %>% 
  chisq_test()
# A tibble: 1 × 6
      n statistic     p    df method          p.signif
* <dbl>     <dbl> <dbl> <int> <chr>           <chr>   
1  5888      3.53 0.473     4 Chi-square test ns      

Có rất nhiều hàm và kiểm định thống kê có thể được thực hiện bằng các hàm trong package rstatix. Đọc các tài liệu về rstatix online ở đây hoặc gõ ?rstatix.

18.4 gtsummary package

Sử dụng package gtsummary nếu bạn đang muốn thêm kết quả của một kiểm định thống kê vào một bảng đẹp được tạo ra bằng package này (như đã được mô tả trong phần gtsummary của chương Bảng mô tả).

Khi thực hiện các kiểm định so sánh bằng hàm tbl_summary, dùng thêm hàm add_p để đưa cột giá trị p và kiểm định được sử dụng vào bảng. Có thể xuất nhiều giá trị p mà được hiệu chỉnh cho nhiều kiểm định bằng cách dùng thêm hàm add_q. Gõ lệnh ?tbl_summary để biết thêm chi tiết.

Kiểm định Chi bình phương

Được sử dụng để so sánh các tỷ lệ của một biến phân loại trong hai nhóm. Kiểm định thống kê mặc định cho biến phân loại trong hàm add_p() là kiểm định Chi bình phương về tính độc lập với hiệu chỉnh liên tục, nhưng nếu có bất kỳ giá trị kỳ vọng nào nhỏ hơn 5 thì kiểm định chính xác của Fisher sẽ được sử dụng.

linelist %>% 
  select(gender, outcome) %>%    # keep variables of interest
  tbl_summary(by = outcome) %>%  # produce summary table and specify grouping variable
  add_p()                        # specify what test to perform
1323 observations missing `outcome` have been removed. To include these observations, use `forcats::fct_na_value_to_level()` on `outcome` column before passing to `tbl_summary()`.

Characteristic

Death, N = 2,582

1

Recover, N = 1,983

1

p-value

2
gender

>0.9
    f 1,227 (50%) 953 (50%)
    m 1,228 (50%) 950 (50%)
    Unknown 127 80
1

n (%)

2

Pearson’s Chi-squared test

Kiểm định t

Được sử dụng để so sánh sự khác biệt về trung bình của một biến trung bình trong hai nhóm. Ví dụ như so sánh tuổi trung bình với kết cục của bệnh nhân.

linelist %>% 
  select(age_years, outcome) %>%             # keep variables of interest
  tbl_summary(                               # produce summary table
    statistic = age_years ~ "{mean} ({sd})", # specify what statistics to show
    by = outcome) %>%                        # specify the grouping variable
  add_p(age_years ~ "t.test")                # specify what tests to perform
1323 observations missing `outcome` have been removed. To include these observations, use `forcats::fct_na_value_to_level()` on `outcome` column before passing to `tbl_summary()`.

Characteristic

Death, N = 2,582

1

Recover, N = 1,983

1

p-value

2
age_years 16 (12) 16 (13) 0.6
    Unknown 32 28
1

Mean (SD)

2

Welch Two Sample t-test

Kiểm định tổng thứ hạng Wilcoxon

Được dùng để so sánh sự phân bố của một biến liên tục trong hai nhóm. Kiểm định mặc định là kiểm định tổng thứ hang Wilcoxon và trung vị (khoảng tứ phân vị IQR) khi so sánh hai nhóm. Tuy nhiên, đối với số liệu không có phân bố chuẩn hoặc so sánh nhiều nhóm, kiểm định Kruskal-wallis là kiểm định thích hợp hơn.

linelist %>% 
  select(age_years, outcome) %>%                       # keep variables of interest
  tbl_summary(                                         # produce summary table
    statistic = age_years ~ "{median} ({p25}, {p75})", # specify what statistic to show (this is default so could remove)
    by = outcome) %>%                                  # specify the grouping variable
  add_p(age_years ~ "wilcox.test")                     # specify what test to perform (default so could leave brackets empty)
1323 observations missing `outcome` have been removed. To include these observations, use `forcats::fct_na_value_to_level()` on `outcome` column before passing to `tbl_summary()`.

Characteristic

Death, N = 2,582

1

Recover, N = 1,983

1

p-value

2
age_years 13 (6, 23) 13 (6, 23) 0.8
    Unknown 32 28
1

Median (IQR)

2

Wilcoxon rank sum test

Kiểm định Kruskal-wallis

Được sử dụng để so sánh sự phân bố của một biến liên tục trong hai hay nhiều nhóm, bất kể số liệu có phân bố chuẩn hay không.

linelist %>% 
  select(age_years, outcome) %>%                       # keep variables of interest
  tbl_summary(                                         # produce summary table
    statistic = age_years ~ "{median} ({p25}, {p75})", # specify what statistic to show (default, so could remove)
    by = outcome) %>%                                  # specify the grouping variable
  add_p(age_years ~ "kruskal.test")                    # specify what test to perform
1323 observations missing `outcome` have been removed. To include these observations, use `forcats::fct_na_value_to_level()` on `outcome` column before passing to `tbl_summary()`.

Characteristic

Death, N = 2,582

1

Recover, N = 1,983

1

p-value

2
age_years 13 (6, 23) 13 (6, 23) 0.8
    Unknown 32 28
1

Median (IQR)

2

Kruskal-Wallis rank sum test

18.5 Tương quan

Mối tương quan giữa các biến định lượng có thể được kiển bằng cách sử dụng lệnh corrr từ package tidyverse. Lệnh này cũng cho phép tính các hệ số tương quan bằng phương pháp Pearson, Kendall hoặc Spearman. Gói lệnh này tạo ra một bảng kết quả và cũng có chức năng tự động vẽ các giá trị.

correlation_tab <- linelist %>% 
  select(generation, age, ct_blood, days_onset_hosp, wt_kg, ht_cm) %>%   # keep numeric variables of interest
  correlate()      # create correlation table (using default pearson)

correlation_tab    # print
# A tibble: 6 × 7
  term            generation      age ct_blood days_onset_hosp    wt_kg    ht_cm
  <chr>                <dbl>    <dbl>    <dbl>           <dbl>    <dbl>    <dbl>
1 generation        NA       -2.22e-2  0.179         -0.288    -0.0302  -0.00942
2 age               -0.0222  NA        0.00849       -0.000635  0.833    0.877  
3 ct_blood           0.179    8.49e-3 NA             -0.600    -0.00636  0.0181 
4 days_onset_hosp   -0.288   -6.35e-4 -0.600         NA         0.0153  -0.00953
5 wt_kg             -0.0302   8.33e-1 -0.00636        0.0153   NA        0.884  
6 ht_cm             -0.00942  8.77e-1  0.0181        -0.00953   0.884   NA      
## remove duplicate entries (the table above is mirrored) 
correlation_tab <- correlation_tab %>% 
  shave()

## view correlation table 
correlation_tab
# A tibble: 6 × 7
  term            generation       age ct_blood days_onset_hosp  wt_kg ht_cm
  <chr>                <dbl>     <dbl>    <dbl>           <dbl>  <dbl> <dbl>
1 generation        NA       NA        NA              NA       NA        NA
2 age               -0.0222  NA        NA              NA       NA        NA
3 ct_blood           0.179    0.00849  NA              NA       NA        NA
4 days_onset_hosp   -0.288   -0.000635 -0.600          NA       NA        NA
5 wt_kg             -0.0302   0.833    -0.00636         0.0153  NA        NA
6 ht_cm             -0.00942  0.877     0.0181         -0.00953  0.884    NA
## plot correlations 
rplot(correlation_tab)

18.6 Nguồn

Phần lớn thông tin trong phần này được phỏng theo các nguồn sau:

gtsummary dplyr corrr sthda correlation