students <- read.csv("students.csv")
df <- students[ , -c(1,3,8,10,12)]
df <- df |> subset(graduated == 1)
row.names(df) <- NULL
# head(df)Kiểm định thống kê so sánh giữa các nhóm
Tài liệu biên soạn trong chuyên đề R Statistics
1 Import dataset
Sử dụng dataset students.csv, chọn những bạn đã tốt nghiệp (graduated == 1).
1.1 Chia group theo độ tuổi
library(dplyr)
df |> dplyr:::mutate(age_group = case_when(
age >= 30 ~ ">= 30",
age < 30 ~ "< 30",
.default = NA
)) -> df
df$age_group <- factor(df$age_group)1.2 Subset 50 bạn trong nhóm tuổi
id_upper_30 <- which(df$age_group == ">= 30")
set.seed(1)
sample(x = id_upper_30,
size = 50,
replace = F) -> id_upper_30_ok
id_upper_30_ok [1] 1486 792 10 700 888 295 1264 1009 470 1193 1573 113 724 1548 1661 904 687 712 866 781 962 434 565 1388 1250 1522 1238 500 91 816 1628
[32] 550 1237 1705 1430 1077 529 658 1637 46 382 492 935 1401 1691 969 1298 1625 180 373id_under_30 <- which(df$age_group == "< 30")
set.seed(1)
sample(x = id_under_30,
size = 50,
replace = F) -> id_under_30_ok
id_under_30_ok [1] 1068 714 137 979 1609 493 313 283 1270 1395 628 1362 1593 344 1697 38 1182 766 923 508 919 401 1708 1026 842 895 980 340 1740 583 1514
[32] 119 1496 558 1606 585 934 358 612 1200 41 564 394 260 1281 397 40 455 851 1481df_sample <- df[ c(id_upper_30_ok, id_under_30_ok) , ]
row.names(df_sample) <- NULL# aggregate(salary ~ age_group,
# data = df_sample,
# FUN = mean)# điều chỉnh dữ liệu cho nhóm upper luôn cao hơn under
set.seed(1)
df_sample$salary[which(df_sample$age_group == ">= 30")] <- df_sample$salary[which(df_sample$age_group == ">= 30")] + 1500df_sample <- df_sample[ , -11]
df_sample$ID <- c(1:length(id_upper_30_ok),
1:length(id_under_30_ok))
dim(df_sample)[1] 100 132 Independent two-sample | So sánh giữa hai nhóm độc lập
2.1 Parametric independent two-sample t-test
2.1.1 Check normality | Kiểm tra giả thuyết phân bố chuẩn
H0: normality (p-value >= 0.05)
H1: NOT normality (p-value < 0.05)
upper_30 <- df_sample$salary[which(df_sample$age_group == ">= 30")]
under_30 <- df_sample$salary[which(df_sample$age_group == "< 30")]library(car)
car:::qqPlot(upper_30)
[1] 38 23car:::qqPlot(under_30)
[1] 33 35shapiro.test(upper_30)
Shapiro-Wilk normality test
data: upper_30
W = 0.96657, p-value = 0.1673ks.test(x = jitter(upper_30),
y = "pnorm",
mean = mean(upper_30),
sd = sd(upper_30))
Exact one-sample Kolmogorov-Smirnov test
data: jitter(upper_30)
D = 0.1036, p-value = 0.6192
alternative hypothesis: two-sidedlibrary(nortest)
nortest:::lillie.test(upper_30)
Lilliefors (Kolmogorov-Smirnov) normality test
data: upper_30
D = 0.10354, p-value = 0.1987nortest:::ad.test(upper_30)
Anderson-Darling normality test
data: upper_30
A = 0.80719, p-value = 0.03421shapiro.test(under_30)
Shapiro-Wilk normality test
data: under_30
W = 0.97963, p-value = 0.5366ks.test(x = jitter(under_30),
y = "pnorm",
mean = mean(under_30),
sd = sd(under_30))
Exact one-sample Kolmogorov-Smirnov test
data: jitter(under_30)
D = 0.073134, p-value = 0.934
alternative hypothesis: two-sidedlibrary(nortest)
nortest:::lillie.test(under_30)
Lilliefors (Kolmogorov-Smirnov) normality test
data: under_30
D = 0.073247, p-value = 0.7217nortest:::ad.test(under_30)
Anderson-Darling normality test
data: under_30
A = 0.37747, p-value = 0.39612.1.2 Check variance homogeneity | Kiểm tra giả thuyết phương sai đồng nhất
H0: phương sai giữa các nhóm bằng nhau
H1: có ít nhất một nhóm có phương sai khác biệt
# kiểm định giữa 2 nhóm
# phổ biến nhất, không yêu cầu 2 vector bằng nhau
stats:::var.test(upper_30,
under_30)
F test to compare two variances
data: upper_30 and under_30
F = 1.0588, num df = 49, denom df = 49, p-value = 0.8424
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.6008322 1.8657694
sample estimates:
ratio of variances
1.05878 stats:::var.test(salary ~ age_group,
data = df_sample)
F test to compare two variances
data: salary by age_group
F = 0.94448, num df = 49, denom df = 49, p-value = 0.8424
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.5359719 1.6643581
sample estimates:
ratio of variances
0.9444835 # kiểm định giữa 2 hay nhiều nhóm
# mạnh khi dữ liệu có phân phối chuẩn
# nhạy với vi phạm giả định chuẩn
stats:::bartlett.test(df_sample$salary, # vector numeric
df_sample$age_group) # vector categorical
Bartlett test of homogeneity of variances
data: df_sample$salary and df_sample$age_group
Bartlett's K-squared = 0.039555, df = 1, p-value = 0.8424# kiểm định giữa 2 hay nhiều nhóm
# phù hợp khi dữ liệu chuẩn và không chuẩn
library(car)
# dữ liệu phân bố chuẩn
car::leveneTest(df_sample$salary,
df_sample$age_group,
center = "mean")Levene's Test for Homogeneity of Variance (center = "mean")
Df F value Pr(>F)
group 1 0.7997 0.3734
98 # dữ liệu lệch chuẩn Brown–Forsythe test (biến thể của Levene's test)
car::leveneTest(df_sample$salary,
df_sample$age_group,
center = "median")Levene's Test for Homogeneity of Variance (center = "median")
Df F value Pr(>F)
group 1 0.5789 0.4486
98 # kiểm định giữa 2 hay nhiều nhóm
# khi dữ liệu không có phân bố chuẩn
# Performs a Fligner-Killeen (median) test of the null that the variances in each of the groups (samples) are the same.
stats::fligner.test(df_sample$salary,
df_sample$major)
Fligner-Killeen test of homogeneity of variances
data: df_sample$salary and df_sample$major
Fligner-Killeen:med chi-squared = 16.767, df = 5, p-value = 0.0049642.1.3 Vẽ đồ thị boxplot
boxplot(upper_30, under_30,
las = 1,
names = c("upper_30", "under_30"),
main = "Salary | upper_30 vs under_30")
2.1.4 Kiểm định two-sample t-test
2.1.4.1 Dữ liệu phân bố chuẩn và phương sai đồng nhất
NORMALITY | VARIANCE HOMOGENEITY (EQUAL VARIANCE)
# TWO-TAILS
# H0: trung binh upper_30 = trung binh under_30
# H1: trung binh upper_30 ≠ trung binh under_30
t.test(upper_30,
under_30,
var.equal = TRUE,
paired = FALSE,
alternative = "two.sided")
Two Sample t-test
data: upper_30 and under_30
t = 2.2423, df = 98, p-value = 0.0272
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
538.9874 8836.1662
sample estimates:
mean of x mean of y
45855.14 41167.57 # LEFT-TAIL
# H0: trung binh upper_30 ≥ trung binh under_30
# H1: trung binh upper_30 < trung binh under_30
t.test(upper_30,
under_30,
var.equal = TRUE,
paired = FALSE,
alternative = "less")
Two Sample t-test
data: upper_30 and under_30
t = 2.2423, df = 98, p-value = 0.9864
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf 8159.009
sample estimates:
mean of x mean of y
45855.14 41167.57 # RIGHT-TAIL
# H0: trung binh upper_30 ≤ trung binh under_30
# H1: trung binh upper_30 > trung binh under_30
t.test(upper_30,
under_30,
var.equal = TRUE,
paired = FALSE,
alternative = "greater")
Two Sample t-test
data: upper_30 and under_30
t = 2.2423, df = 98, p-value = 0.0136
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
1216.144 Inf
sample estimates:
mean of x mean of y
45855.14 41167.57 2.1.4.2 Dữ liệu phân bố chuẩn và phương sai không đồng nhất
NORMALITY and VARIANCE HETEROSKEDASTICITY (UNEQUAL VARIANCE)
Sử dụng t-test cải tiến | Welch Two Sample t-test
t.test(upper_30,
under_30,
var.equal = FALSE,
paired = FALSE,
alternative = "two.sided")
Welch Two Sample t-test
data: upper_30 and under_30
t = 2.2423, df = 97.92, p-value = 0.0272
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
538.9451 8836.2085
sample estimates:
mean of x mean of y
45855.14 41167.57 2.1.4.3 Effect size
Compute the effect size for t-test. T-test conventional effect sizes, proposed by Cohen, are: 0.2 (small effect), 0.5 (moderate effect) and 0.8 (large effect).
library(rstatix)
rstatix:::cohens_d(salary ~ age_group,
data = df_sample,
var.equal = FALSE,
paired = FALSE)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 salary < 30 >= 30 -0.448 50 50 small library(effectsize)
effectsize:::cohens_d(
upper_30,
under_30,
paired = F
)Cohen's d | 95% CI
------------------------
0.45 | [0.05, 0.84]
- Estimated using pooled SD.2.2 Nonparametric
2.2.1 Independent two-sample Mann–Whitney U test
Mann–Whitney U test (aka Wilcoxon Rank Sum test) compares two independent groups. Dành cho dữ liệu không đạt phân phối chuẩn (gần chuẩn).
Nonparametric không yêu cầu giả định NORMALITY và VARIANCE HOMOGENEITY.
wilcox.test(upper_30,
under_30,
paired = FALSE,
alternative = "two.sided")
Wilcoxon rank sum test with continuity correction
data: upper_30 and under_30
W = 1563, p-value = 0.03122
alternative hypothesis: true location shift is not equal to 0library(rstatix)
# nếu chỉ có so sánh giữa 2 nhóm thì không cần adjust p-value
rstatix:::wilcox_test(salary ~ age_group,
data = df_sample,
paired = FALSE,
p.adjust.method = "none",
alternative = "two.sided") # A tibble: 1 × 7
.y. group1 group2 n1 n2 statistic p
* <chr> <chr> <chr> <int> <int> <dbl> <dbl>
1 salary < 30 >= 30 50 50 937 0.03122.2.2 Effect size
The r value varies from 0 to close to 1. The interpretation values for r commonly in published litterature and on the internet are: 0.10 - < 0.3 (small effect), 0.30 - < 0.5 (moderate effect) and >= 0.5 (large effect).
rstatix:::wilcox_effsize(salary ~ age_group,
data = df_sample,
paired = FALSE)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 salary < 30 >= 30 0.216 50 50 small 2.2.3 Independent two-sample Mood’s Median Test
Performs a Mood’s median test to compare medians of independent samples
# https://rcompanion.org/handbook/F_09.html
library(RVAideMemoire)
RVAideMemoire:::mood.medtest(salary ~ age_group,
data = df_sample,
exact = FALSE)
Mood's median test
data: salary by age_group
X-squared = 1.96, df = 1, p-value = 0.1615Testing the equality of the distributions of a numeric response variable in two or more independent groups against shift alternatives.
library(coin)
coin:::median_test(salary ~ age_group,
data = df_sample)
Asymptotic Two-Sample Brown-Mood Median Test
data: salary by age_group (< 30, >= 30)
Z = -1.592, p-value = 0.1114
alternative hypothesis: true mu is not equal to 02.3 Áp dụng package
library(ggstatsplot)
library(tidyverse)
set.seed(20)
ggstatsplot:::ggbetweenstats(data = df_sample,
x = age_group,
y = salary,
type = "parametric",
var.equal = TRUE,
plot.type = "box",
pairwise.comparisons = T,
pairwise.display = "all",
p.adjust.method = "none",
centrality.plotting = T,
centrality.path = T,
digits = 4L,
centrality.point.args = list(size = 4, color = "darkred"),
centrality.label.args = list(size = 4,
nudge_x = 0.4,
segment.linetype = 4),
ggtheme = theme_classic(base_size = 12),
bf.message = T,
results.subtitle = T,
caption = "",
ggsignif.args = list(size = 0.5,
vjust = 0,
step_increase = 0.1,
textsize = 5),
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6),
alpha = 0.5,
size = 6,
stroke = 0,
na.rm = TRUE)
) 
library(dabestr)
two_groups_unpaired <- dabestr:::load(df_sample,
x = age_group, y = salary,
paired = NULL,
idx = c("< 30", ">= 30")
)
two_groups_unpaired.mean_diff <- dabestr::cohens_d(two_groups_unpaired)
print(two_groups_unpaired.mean_diff)DABESTR v2025.3.15
==================
Good evening!
The current time is 23:02 PM on Thursday October 29, 2026.
The character(0) Cohen's d between >= 30 and < 30 is 0.448 [95%CI 0.054, 0.851].
The p-value of the two-sided permutation t-test is 0.0272, calculated for legacy purposes only.
5000 bootstrap samples were taken; the confidence interval is bias-corrected and accelerated.
Any p-value reported is the probability of observing the effect size (or greater),
assuming the null hypothesis of zero difference is true.
For each p-value, 5000 reshuffles of the control and test labels were performed.dabest_plot(two_groups_unpaired.mean_diff)
https://rpubs.com/lengockhanhi/454078
3 Dependent two sample
3.1 Parametric paired t-test
t.test(df_sample$score1,
df_sample$score2,
var.equal = TRUE,
paired = TRUE,
alternative = "two.sided")
Paired t-test
data: df_sample$score1 and df_sample$score2
t = -1.2863, df = 99, p-value = 0.2014
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-2.2629386 0.4829386
sample estimates:
mean difference
-0.89 t.test(df_sample$score1 - df_sample$score2,
mu = 0,
alternative = "two.sided")
One Sample t-test
data: df_sample$score1 - df_sample$score2
t = -1.2863, df = 99, p-value = 0.2014
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-2.2629386 0.4829386
sample estimates:
mean of x
-0.89 # head(df_sample)
df_score <- df_sample[ , c("stud.id", "score1", "score2")]
df_score stud.id score1 score2
1 127329 71 73
2 699675 69 59
3 373120 56 52
4 645027 60 48
5 887121 81 73
6 263680 74 64
7 644535 90 85
8 920111 85 88
9 534172 68 73
10 152383 54 45
11 165588 67 77
12 358961 93 90
13 854188 71 63
14 355955 50 59
15 684695 55 63
16 526849 59 53
17 731560 74 69
18 728977 76 83
19 244176 87 93
20 337141 72 73
21 714747 70 79
22 245697 70 77
23 606531 49 58
24 149335 74 70
25 281483 76 62
26 195794 47 55
27 488797 83 78
28 451408 73 74
29 729441 60 57
30 587073 58 60
31 195357 40 48
32 863615 72 73
33 219503 51 57
34 551154 78 83
35 672951 81 90
36 969474 74 81
37 831973 74 80
38 364152 85 87
39 287337 92 89
40 192380 90 92
41 970464 86 90
42 383269 56 57
43 597038 86 80
44 141917 64 49
45 885855 74 68
46 906012 59 59
47 553535 51 51
48 456353 65 63
49 654954 83 89
50 988603 69 68
51 253803 67 72
52 595839 66 60
53 711793 63 54
54 429789 52 44
55 130608 48 53
56 461293 67 68
57 448938 90 86
58 809408 70 69
59 529799 84 86
60 786487 88 88
61 918215 90 94
62 447520 72 78
63 859275 76 78
64 861582 88 86
65 394142 79 76
66 595769 78 77
67 846250 67 79
68 747054 72 86
69 749611 58 51
70 748309 85 90
71 139052 42 40
72 921795 46 62
73 869910 61 54
74 678552 61 69
75 520250 61 80
76 819078 36 51
77 371542 67 74
78 122995 78 76
79 565624 66 80
80 350085 59 52
81 648615 85 86
82 621166 73 69
83 760133 81 74
84 211707 74 76
85 320119 61 56
86 249763 73 66
87 181738 75 65
88 707668 75 74
89 111998 64 72
90 869321 73 75
91 612668 55 64
92 685552 65 70
93 726881 56 53
94 288990 85 87
95 681749 67 77
96 923689 59 61
97 390559 78 83
98 954089 57 57
99 901358 75 79
100 160166 81 76###
df_score1 <- df_sample[ , c("stud.id", "score1")]
names(df_score1) <- c("id", "score")
df_score1$group <- "hk1"
# df_score1
###
df_score2 <- df_sample[ , c("stud.id", "score2")]
names(df_score2) <- c("id", "score")
df_score2$group <- "hk2"
# df_score2
df_long <- rbind(df_score1,
df_score2)
df_long id score group
1 127329 71 hk1
2 699675 69 hk1
3 373120 56 hk1
4 645027 60 hk1
5 887121 81 hk1
6 263680 74 hk1
7 644535 90 hk1
8 920111 85 hk1
9 534172 68 hk1
10 152383 54 hk1
11 165588 67 hk1
12 358961 93 hk1
13 854188 71 hk1
14 355955 50 hk1
15 684695 55 hk1
16 526849 59 hk1
17 731560 74 hk1
18 728977 76 hk1
19 244176 87 hk1
20 337141 72 hk1
21 714747 70 hk1
22 245697 70 hk1
23 606531 49 hk1
24 149335 74 hk1
25 281483 76 hk1
26 195794 47 hk1
27 488797 83 hk1
28 451408 73 hk1
29 729441 60 hk1
30 587073 58 hk1
31 195357 40 hk1
32 863615 72 hk1
33 219503 51 hk1
34 551154 78 hk1
35 672951 81 hk1
36 969474 74 hk1
37 831973 74 hk1
38 364152 85 hk1
39 287337 92 hk1
40 192380 90 hk1
41 970464 86 hk1
42 383269 56 hk1
43 597038 86 hk1
44 141917 64 hk1
45 885855 74 hk1
46 906012 59 hk1
47 553535 51 hk1
48 456353 65 hk1
49 654954 83 hk1
50 988603 69 hk1
51 253803 67 hk1
52 595839 66 hk1
53 711793 63 hk1
54 429789 52 hk1
55 130608 48 hk1
56 461293 67 hk1
57 448938 90 hk1
58 809408 70 hk1
59 529799 84 hk1
60 786487 88 hk1
61 918215 90 hk1
62 447520 72 hk1
63 859275 76 hk1
64 861582 88 hk1
65 394142 79 hk1
66 595769 78 hk1
67 846250 67 hk1
68 747054 72 hk1
69 749611 58 hk1
70 748309 85 hk1
71 139052 42 hk1
72 921795 46 hk1
73 869910 61 hk1
74 678552 61 hk1
75 520250 61 hk1
76 819078 36 hk1
77 371542 67 hk1
78 122995 78 hk1
79 565624 66 hk1
80 350085 59 hk1
81 648615 85 hk1
82 621166 73 hk1
83 760133 81 hk1
84 211707 74 hk1
85 320119 61 hk1
86 249763 73 hk1
87 181738 75 hk1
88 707668 75 hk1
89 111998 64 hk1
90 869321 73 hk1
91 612668 55 hk1
92 685552 65 hk1
93 726881 56 hk1
94 288990 85 hk1
95 681749 67 hk1
96 923689 59 hk1
97 390559 78 hk1
98 954089 57 hk1
99 901358 75 hk1
100 160166 81 hk1
101 127329 73 hk2
102 699675 59 hk2
103 373120 52 hk2
104 645027 48 hk2
105 887121 73 hk2
106 263680 64 hk2
107 644535 85 hk2
108 920111 88 hk2
109 534172 73 hk2
110 152383 45 hk2
111 165588 77 hk2
112 358961 90 hk2
113 854188 63 hk2
114 355955 59 hk2
115 684695 63 hk2
116 526849 53 hk2
117 731560 69 hk2
118 728977 83 hk2
119 244176 93 hk2
120 337141 73 hk2
121 714747 79 hk2
122 245697 77 hk2
123 606531 58 hk2
124 149335 70 hk2
125 281483 62 hk2
126 195794 55 hk2
127 488797 78 hk2
128 451408 74 hk2
129 729441 57 hk2
130 587073 60 hk2
131 195357 48 hk2
132 863615 73 hk2
133 219503 57 hk2
134 551154 83 hk2
135 672951 90 hk2
136 969474 81 hk2
137 831973 80 hk2
138 364152 87 hk2
139 287337 89 hk2
140 192380 92 hk2
141 970464 90 hk2
142 383269 57 hk2
143 597038 80 hk2
144 141917 49 hk2
145 885855 68 hk2
146 906012 59 hk2
147 553535 51 hk2
148 456353 63 hk2
149 654954 89 hk2
150 988603 68 hk2
151 253803 72 hk2
152 595839 60 hk2
153 711793 54 hk2
154 429789 44 hk2
155 130608 53 hk2
156 461293 68 hk2
157 448938 86 hk2
158 809408 69 hk2
159 529799 86 hk2
160 786487 88 hk2
161 918215 94 hk2
162 447520 78 hk2
163 859275 78 hk2
164 861582 86 hk2
165 394142 76 hk2
166 595769 77 hk2
167 846250 79 hk2
168 747054 86 hk2
169 749611 51 hk2
170 748309 90 hk2
171 139052 40 hk2
172 921795 62 hk2
173 869910 54 hk2
174 678552 69 hk2
175 520250 80 hk2
176 819078 51 hk2
177 371542 74 hk2
178 122995 76 hk2
179 565624 80 hk2
180 350085 52 hk2
181 648615 86 hk2
182 621166 69 hk2
183 760133 74 hk2
184 211707 76 hk2
185 320119 56 hk2
186 249763 66 hk2
187 181738 65 hk2
188 707668 74 hk2
189 111998 72 hk2
190 869321 75 hk2
191 612668 64 hk2
192 685552 70 hk2
193 726881 53 hk2
194 288990 87 hk2
195 681749 77 hk2
196 923689 61 hk2
197 390559 83 hk2
198 954089 57 hk2
199 901358 79 hk2
200 160166 76 hk2p1 <- ggstatsplot:::ggwithinstats(data = df_long,
x = group,
y = score,
type = "parametric",
var.equal = TRUE,
plot.type = "box",
pairwise.comparisons = T,
pairwise.display = "all",
p.adjust.method = "none",
centrality.plotting = T,
centrality.path = T,
digits = 4L,
ggtheme = theme_classic(base_size = 12),
bf.message = T,
results.subtitle = T,
caption = "",
centrality.point.args = list(size = 4, color = "darkred"),
centrality.label.args = list(size = 4,
nudge_x = 0.4,
segment.linetype = 4),
ggsignif.args = list(size = 0.5,
vjust = 0,
step_increase = 0.1,
textsize = 5),
point.path.args = list(alpha = 0.5,
lwd = 0.8,
linetype = "solid",
color = "blue"),
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0, jitter.width = 0, jitter.height = 0),
alpha = 0.5,
size = 6,
stroke = 0,
na.rm = TRUE)
)
ggsave(filename = "p1_paired.png",
width = 8,
height = 6,
plot = p1,
dpi = 300,
units = "in")
library(dabestr)
two_groups_paired <- dabestr:::load(df_long,
x = group,
y = score,
paired = "sequential",
idx = c("hk1", "hk2"),
id_col = id
)
two_groups_paired.mean_diff <- dabestr:::mean_diff(two_groups_paired)
print(two_groups_paired.mean_diff)DABESTR v2025.3.15
==================
Good evening!
The current time is 23:02 PM on Thursday October 29, 2026.
The paired mean difference between hk2 and hk1 is 0.89 [95%CI -2.85, 4.404].
The p-value of the two-sided permutation t-test is 0.5916, calculated for legacy purposes only.
5000 bootstrap samples were taken; the confidence interval is bias-corrected and accelerated.
Any p-value reported is the probability of observing the effect size (or greater),
assuming the null hypothesis of zero difference is true.
For each p-value, 5000 reshuffles of the control and test labels were performed.dabest_plot(two_groups_paired.mean_diff)
library(rstatix)
rstatix:::cohens_d(salary ~ age_group,
data = df_sample,
var.equal = T,
paired = T)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 salary < 30 >= 30 -0.310 50 50 small library(effectsize)
effectsize:::cohens_d(
upper_30,
under_30,
paired = T
)Cohen's d | 95% CI
------------------------
0.31 | [0.02, 0.59]3.2 Nonparametric | Paired
3.2.1 Wilcoxon Signed Rank test
wilcox.test(upper_30,
under_30,
paired = T,
alternative = "two.sided")
Wilcoxon signed rank test with continuity correction
data: upper_30 and under_30
V = 858, p-value = 0.03369
alternative hypothesis: true location shift is not equal to 03.2.2 Sign test for paired samples
Khi phân phối lệch chuẩn quá nhiều thì ta dùng sign test.
library(BSDA)
BSDA:::SIGN.test(upper_30,
under_30,
md = 0,
alternative = "two.sided",
conf.level = 0.95)
Dependent-samples Sign-Test
data: upper_30 and under_30
S = 33, p-value = 0.03284
alternative hypothesis: true median difference is not equal to 0
95 percent confidence interval:
803.5553 9798.7770
sample estimates:
median of x-y
4115.385
Achieved and Interpolated Confidence Intervals:
Conf.Level L.E.pt U.E.pt
Lower Achieved CI 0.9351 955.7794 9699.373
Interpolated CI 0.9500 803.5553 9798.777
Upper Achieved CI 0.9672 628.3350 9913.197library(DescTools)
DescTools:::SignTest(upper_30,
under_30,
md = 0,
alternative = "two.sided",
conf.level = 0.95)
Dependent-samples Sign-Test
data: upper_30 and under_30
S = 33, number of differences = 50, p-value = 0.03284
alternative hypothesis: true median difference is not equal to 0
96.7 percent confidence interval:
628.335 9913.197
sample estimates:
median of the differences
4115.385 4 t-test cho nhiều nhóm
library(ggplot2)
library(ggsignif)
options(scipen = 999)
dp <- df[ , c("gender",
"major",
"salary")]
dim(dp)[1] 1753 3dp <- na.omit(dp)
dim(dp)[1] 1753 3p1 <- ggplot(data = dp,
mapping = aes(x = gender,
y = salary,
fill = gender)) +
# geom_boxplot() +
geom_violin(alpha = 0.1,
show.legend = F) +
geom_point(position = position_jitter(width = 0.1, height = 0),
mapping = aes(color = gender),
alpha = 0.5) +
stat_summary(fun.data = mean_cl_boot,
geom ='crossbar',
width = 0.6,
alpha = 0.5,
aes(fill = gender)) +
ggsignif:::geom_signif(comparisons = list(c("Male", "Female")),
map_signif_level = function(p_ok) {
options(scipen = 1)
options(digits = 2)
capture.output(cat(p_ok)) -> p_ok_final
paste0("p", " = ", p_ok_final)
},
textsize = 4,
vjust = -0.5,
step_increase = 0.6,
test = "t.test",
test.args = list(paired = F)
) +
facet_grid(~ major) +
scale_y_continuous(limits = c(10000,100000),
breaks = seq(from = 10000,
to = 100000,
by = 10000)) +
theme_bw(base_size = 13) +
labs(x = "")
# p1
ggsave(filename = "p1.png",
plot = p1,
width = 16,
height = 4,
units = "in",
dpi = 300)
5 ANOVA 1 yếu tố | Parametric | Không bắt cặp
5.1 Check assumption
5.1.1 Normality
library(ggplot2)
library(ggpubr)
car:::qqPlot(PlantGrowth$weight[1:10])
[1] 4 1ggpubr:::ggqqplot(data = PlantGrowth,
x = "weight",
facet.by = "group")
ctrl <- PlantGrowth$weight[1:10]
shapiro.test(ctrl)
Shapiro-Wilk normality test
data: ctrl
W = 1, p-value = 0.7ks.test(x = jitter(ctrl),
y = "pnorm",
mean = mean(ctrl),
sd = sd(ctrl))
Exact one-sample Kolmogorov-Smirnov test
data: jitter(ctrl)
D = 0.2, p-value = 0.9
alternative hypothesis: two-sidedlibrary(nortest)
nortest:::lillie.test(ctrl)
Lilliefors (Kolmogorov-Smirnov) normality test
data: ctrl
D = 0.2, p-value = 0.5nortest:::ad.test(ctrl)
Anderson-Darling normality test
data: ctrl
A = 0.3, p-value = 0.6trt1 <- PlantGrowth$weight[11:20]
shapiro.test(trt1)
Shapiro-Wilk normality test
data: trt1
W = 0.9, p-value = 0.5ks.test(x = jitter(trt1),
y = "pnorm",
mean = mean(trt1),
sd = sd(trt1))
Exact one-sample Kolmogorov-Smirnov test
data: jitter(trt1)
D = 0.2, p-value = 0.8
alternative hypothesis: two-sidedlibrary(nortest)
nortest:::lillie.test(trt1)
Lilliefors (Kolmogorov-Smirnov) normality test
data: trt1
D = 0.2, p-value = 0.4nortest:::ad.test(trt1)
Anderson-Darling normality test
data: trt1
A = 0.3, p-value = 0.4library(tidyr)
library(dplyr)
# reshape tu long qua wide
PlantGrowth_wide <- PlantGrowth %>%
group_by(group) %>%
mutate(id = row_number()) %>%
ungroup() %>% as.data.frame() %>%
pivot_wider(
names_from = group,
values_from = weight
)
PlantGrowth_wide <- PlantGrowth_wide[, -1]
library(MVN)
MVN:::mvn(PlantGrowth_wide)$multivariate_normality
Test Statistic p.value Method MVN
1 Henze-Zirkler 0.5 0.44 asymptotic ✓ Normal
$univariate_normality
Test Variable Statistic p.value Normality
1 Anderson-Darling ctrl 0.28 0.56 ✓ Normal
2 Anderson-Darling trt1 0.33 0.44 ✓ Normal
3 Anderson-Darling trt2 0.31 0.50 ✓ Normal
$descriptives
Variable n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
1 ctrl 10 5.0 0.58 5.2 4.2 6.1 4.5 5.3 0.27 2.3
2 trt1 10 4.7 0.79 4.5 3.6 6.0 4.2 4.9 0.56 2.3
3 trt2 10 5.5 0.44 5.4 4.9 6.3 5.3 5.7 0.57 2.3
$data
# A tibble: 10 × 3
ctrl trt1 trt2
<dbl> <dbl> <dbl>
1 4.17 4.81 6.31
2 5.58 4.17 5.12
3 5.18 4.41 5.54
4 6.11 3.59 5.5
5 4.5 5.87 5.37
6 4.61 3.83 5.29
7 5.17 6.03 4.92
8 4.53 4.89 6.15
9 5.33 4.32 5.8
10 5.14 4.69 5.26
$subset
NULL
$outlierMethod
[1] "none"
attr(,"class")
[1] "mvn"5.1.2 Variance homogeneity
library(car)
car:::leveneTest(weight ~ group,
data = PlantGrowth)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 1.12 0.34
27 stats:::bartlett.test(weight ~ group,
data = PlantGrowth)
Bartlett test of homogeneity of variances
data: weight by group
Bartlett's K-squared = 3, df = 2, p-value = 0.2library(ggplot2)
ggplot(data = PlantGrowth,
mapping = aes(x = weight,
colour = group)) +
geom_density() +
theme_classic(base_size = 16)
ggplot(data = PlantGrowth,
mapping = aes(x = group,
y = weight,
fill = group)) +
geom_boxplot() +
theme_classic(base_size = 16)
5.2 Run ANOVA one-way
# ANOVA type I | variance equal
fit1 <- aov(weight ~ group,
data = PlantGrowth)
summary(fit1) Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.77 1.883 4.85 0.016 *
Residuals 27 10.49 0.389
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit2 <- stats:::oneway.test(weight ~ group,
data = PlantGrowth,
var.equal = T)
fit2
One-way analysis of means
data: weight and group
F = 5, num df = 2, denom df = 27, p-value = 0.02# ANOVA type I | variance unequal Welch's ANOVA
fit3 <- stats:::oneway.test(weight ~ group,
data = PlantGrowth,
var.equal = F)
fit3
One-way analysis of means (not assuming equal variances)
data: weight and group
F = 5, num df = 2, denom df = 17, p-value = 0.02library(car)
# ANOVA type II
car:::Anova(fit1,
type = "II")Anova Table (Type II tests)
Response: weight
Sum Sq Df F value Pr(>F)
group 3.77 2 4.85 0.016 *
Residuals 10.49 27
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA type III
car:::Anova(fit1, type = "III")Anova Table (Type III tests)
Response: weight
Sum Sq Df F value Pr(>F)
(Intercept) 253.2 1 651.60 <2e-16 ***
group 3.8 2 4.85 0.016 *
Residuals 10.5 27
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.3 Effect size
library(effectsize)
# Eta-squared
effectsize:::eta_squared(fit1)# Effect Size for ANOVA
Parameter | Eta2 | 95% CI
-------------------------------
group | 0.26 | [0.04, 1.00]
- One-sided CIs: upper bound fixed at [1.00].# eta-squared
# Tỷ lệ phương sai của weight được giải thích bởi group.
# 0.01 Nhỏ
# 0.06 Trung bình
# 0.14 Lớn# Omega-squared
# Ước lượng effect size ít thiên lệch hơn trong quần thể.
effectsize:::omega_squared(fit1)# Effect Size for ANOVA
Parameter | Omega2 | 95% CI
---------------------------------
group | 0.20 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].# Cohen's f
# Thước đo effect size dựa trên độ lớn khác biệt giữa các trung bình nhóm.
effectsize:::cohens_f(fit1)# Effect Size for ANOVA
Parameter | Cohen's f | 95% CI
-----------------------------------
group | 0.60 | [0.19, Inf]
- One-sided CIs: upper bound fixed at [Inf].library(rstatix)
# Eta-squared
rstatix:::anova_test(
data = PlantGrowth,
dv = weight,
effect.size = "pes",
between = group)ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 pes
1 group 2 27 4.8 0.016 * 0.265.4 Post-hoc analysis | Multiple comparisons
5.4.1 Tukey
TukeyHSD(fit1) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = PlantGrowth)
$group
diff lwr upr p adj
trt1-ctrl -0.37 -1.06 0.32 0.39
trt2-ctrl 0.49 -0.20 1.19 0.20
trt2-trt1 0.87 0.17 1.56 0.01plot(TukeyHSD(fit1),
las = 1)
library(agricolae)
agricolae:::HSD.test(fit1,
trt = "group",
console = TRUE)
Study: fit1 ~ "group"
HSD Test for weight
Mean Square Error: 0.39
group, means
weight std r se Min Max Q25 Q50 Q75
ctrl 5.0 0.58 10 0.2 4.2 6.1 4.5 5.2 5.3
trt1 4.7 0.79 10 0.2 3.6 6.0 4.2 4.6 4.9
trt2 5.5 0.44 10 0.2 4.9 6.3 5.3 5.4 5.7
Alpha: 0.05 ; DF Error: 27
Critical Value of Studentized Range: 3.5
Minimun Significant Difference: 0.69
Treatments with the same letter are not significantly different.
weight groups
trt2 5.5 a
ctrl 5.0 ab
trt1 4.7 b5.4.2 Bonferroni
• “bonferroni” Kiểm soát family-wise error rate (FWER), khá bảo thủ
• “holm” Kiểm soát FWER, thường ít bảo thủ hơn Bonferroni
• “BH” Kiểm soát false discovery rate (FDR)
pairwise.t.test(
x = PlantGrowth$weight,
g = PlantGrowth$group,
p.adjust.method = "bonferroni",
pool.sd = TRUE
)
Pairwise comparisons using t tests with pooled SD
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.58 -
trt2 0.26 0.01
P value adjustment method: bonferroni pairwise.t.test(
x = PlantGrowth$weight,
g = PlantGrowth$group,
p.adjust.method = "BH",
pool.sd = TRUE
)
Pairwise comparisons using t tests with pooled SD
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.19 -
trt2 0.13 0.01
P value adjustment method: BH pairwise.t.test(
x = PlantGrowth$weight,
g = PlantGrowth$group,
p.adjust.method = "holm",
pool.sd = TRUE
)
Pairwise comparisons using t tests with pooled SD
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.19 -
trt2 0.18 0.01
P value adjustment method: holm pairwise.t.test(
x = PlantGrowth$weight,
g = PlantGrowth$group,
p.adjust.method = "bonferroni",
pool.sd = FALSE
)
Pairwise comparisons using t tests with non-pooled SD
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.75 -
trt2 0.14 0.03
P value adjustment method: bonferroni 5.5 Vẽ đồ thị
library(ggstatsplot)
library(tidyverse)
set.seed(20)
kq <- ggstatsplot:::ggbetweenstats(data = PlantGrowth,
x = group,
y = weight,
var.equal = TRUE,
type = "parametric",
plot.type = "box",
pairwise.comparisons = T,
pairwise.display = "all",
p.adjust.method = "bonferroni",
# p.adjust.method = "none",
centrality.plotting = T,
centrality.path = T,
digits = 4L,
centrality.point.args = list(size = 4, color = "darkred"),
centrality.label.args = list(size = 4,
nudge_x = 0.4,
segment.linetype = 4),
ggtheme = theme_classic(base_size = 12),
bf.message = T,
results.subtitle = T,
caption = "",
ggsignif.args = list(size = 0.5,
vjust = 0,
step_increase = 0.1,
textsize = 5),
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0, jitter.width = 0, jitter.height = 0),
alpha = 0.5,
size = 6,
stroke = 0,
na.rm = TRUE)
) +
coord_cartesian(ylim = c(3,8))
# kq
ggstatsplot:::extract_stats(kq)$subtitle_data
# A tibble: 1 × 14
statistic df df.error p.value method effectsize estimate conf.level conf.low conf.high conf.method conf.distribution n.obs expression
<dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <int> <list>
1 4.85 2 27 0.0159 One-way analysis of means Omega2 0.204 0.95 0.00238 1 ncp F 30 <language>
$caption_data
# A tibble: 6 × 17
term pd prior.distribution prior.location prior.scale bf10 method log_e_bf10 effectsize estimate std.dev conf.level
<chr> <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
1 mu 1 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
2 group-ctrl 0.596 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
3 group-trt1 0.985 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
4 group-trt2 0.989 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
5 sig2 1 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
6 g_group 1 cauchy 0 0.707 3.66 Bayes factors for linear models 1.30 Bayesian R-squared 0.167 0.161 0.95
conf.low conf.high conf.method n.obs expression
<dbl> <dbl> <chr> <int> <list>
1 0 0.381 HDI 30 <language>
2 0 0.381 HDI 30 <language>
3 0 0.381 HDI 30 <language>
4 0 0.381 HDI 30 <language>
5 0 0.381 HDI 30 <language>
6 0 0.381 HDI 30 <language>
$pairwise_comparisons_data
# A tibble: 3 × 6
group1 group2 p.value p.adjust.method test expression
<chr> <chr> <dbl> <chr> <chr> <list>
1 ctrl trt1 0.583 Bonferroni Student's t <language>
2 ctrl trt2 0.263 Bonferroni Student's t <language>
3 trt1 trt2 0.0134 Bonferroni Student's t <language>
$descriptive_data
NULL
$one_sample_data
NULL
$tidy_data
NULL
$glance_data
NULL
attr(,"class")
[1] "ggstatsplot_stats" "list" ggsave(filename = "kq.png",
plot = kq,
width = 8,
height = 6,
units = "in",
dpi = 300)
5.6 Estimated marginal means
Estimated marginal means (EMMs), a.k.a. least-squares means
# https://cran.r-project.org/web/packages/emmeans/vignettes/FAQs.html#what
library(emmeans)
emm <- emmeans:::emmeans(fit1, ~ group)
summary(emm) group emmean SE df lower.CL upper.CL
ctrl 5.0 0.197 27 4.6 5.4
trt1 4.7 0.197 27 4.3 5.1
trt2 5.5 0.197 27 5.1 5.9
Confidence level used: 0.95 graphics:::pairs(emm, adjust = "bonferroni") contrast estimate SE df t.ratio p.value
ctrl - trt1 0.37 0.279 27 1.330 0.5800
ctrl - trt2 -0.49 0.279 27 -1.770 0.2600
trt1 - trt2 -0.86 0.279 27 -3.100 0.0100
P value adjustment: bonferroni method for 3 tests # plot(emm)
# dấu mũi tên 2 chiều là không khác biệt
# hai dấu mũi tên ngược nhau là khác biệt
emmeans:::plot.emmGrid(emm,
comparisons = F)
6 ANOVA 1 yếu tố | Nonparametric | Không bắt cặp
6.1 Check assumption
6.1.1 Normality
library(ggplot2)
library(ggpubr)
df$major <- factor(df$major)
ggpubr:::ggqqplot(data = df,
x = "salary",
facet.by = "major")
major_ok <- unique(df$major)
for(i in 1:length(major_ok)){
check_normality <- df |> subset(major == major_ok[i])
print(major_ok[i])
kq <- nortest:::ad.test(check_normality$salary)
print(kq)
}[1] Economics and Finance
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 0.6, p-value = 0.1
[1] Environmental Sciences
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 0.4, p-value = 0.3
[1] Political Science
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 0.9, p-value = 0.03
[1] Mathematics and Statistics
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 0.8, p-value = 0.05
[1] Biology
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 1, p-value = 0.001
[1] Social Sciences
Levels: Biology Economics and Finance Environmental Sciences Mathematics and Statistics Political Science Social Sciences
Anderson-Darling normality test
data: check_normality$salary
A = 0.3, p-value = 0.6### Lưu ý nếu dữ liệu không đều nhau thì nên thực hiện thủ công trích từng vector kiểm riêng normality
df_wide <- df[ , c("major", "salary")] %>%
group_by(major) %>%
mutate(id = row_number()) %>%
ungroup() %>% as.data.frame() %>%
pivot_wider(
names_from = major,
values_from = salary
)
df_wide <- df_wide[, -1]
library(MVN)
### Vì dữ liệu df_wide bị NA nên khi thực hiện, package sẽ tự động remove NA, làm sai kết quả.
df_wide <- na.omit(df_wide)
MVN:::mvn(df_wide)$multivariate_normality
Test Statistic p.value Method MVN
1 Henze-Zirkler 0.95 0.34 asymptotic ✓ Normal
$univariate_normality
Test Variable Statistic p.value Normality
1 Anderson-Darling Economics and Finance 0.49 0.218 ✓ Normal
2 Anderson-Darling Environmental Sciences 0.14 0.971 ✓ Normal
3 Anderson-Darling Political Science 0.72 0.060 ✓ Normal
4 Anderson-Darling Mathematics and Statistics 0.62 0.102 ✓ Normal
5 Anderson-Darling Biology 0.90 0.021 ✗ Not normal
6 Anderson-Darling Social Sciences 0.30 0.585 ✓ Normal
$descriptives
Variable n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
1 Economics and Finance 206 49760 8423 49349 25351 72787 43373 56022 0.209 2.8
2 Environmental Sciences 206 37229 5681 37132 23870 51332 33221 41368 0.031 2.6
3 Political Science 206 36876 9296 36514 11444 70932 30272 42177 0.480 3.7
4 Mathematics and Statistics 206 48594 8203 47554 23867 75597 43098 53688 0.328 3.3
5 Biology 206 47462 8330 46567 25998 70095 41099 53279 0.279 2.7
6 Social Sciences 206 32044 6205 32032 13300 50705 27613 36021 0.078 3.3
$data
# A tibble: 206 × 6
`Economics and Finance` `Environmental Sciences` `Political Science` `Mathematics and Statistics` Biology `Social Sciences`
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 45254. 40553. 33969. 35599. 49111. 28128.
2 50618. 27007. 37392. 39343. 51223. 23717.
3 42593. 39704. 30087. 36204. 49476. 35882.
4 60407. 32519. 31494. 58081. 46243. 37036.
5 35881. 33069. 39451. 56996. 39588. 39758.
6 51564. 37804. 20511. 41619. 64486. 21722.
7 43380. 35192. 44028. 41373. 45504. 33109.
8 71468. 34191. 37810. 65013. 39369. 31750.
9 51518. 39038. 42065. 64525. 35816. 31485.
10 55664. 43333. 30612. 58671. 37626. 25659.
# ℹ 196 more rows
$subset
NULL
$outlierMethod
[1] "none"
attr(,"class")
[1] "mvn"6.1.2 Variance homogeneity
library(car)
car:::leveneTest(salary ~ major,
data = df)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 16.3 9.6e-16 ***
1747
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1stats:::bartlett.test(salary ~ major,
data = df)
Bartlett test of homogeneity of variances
data: salary by major
Bartlett's K-squared = 108, df = 5, p-value <2e-166.2 Run ANOVA one-way Nonparametric
Assumptions cho parametric bị vi phạm nên ta áp dụng kiểm định nonparametric.
fit1 <- stats:::kruskal.test(salary ~ major,
data = df)
fit1
Kruskal-Wallis rank sum test
data: salary by major
Kruskal-Wallis chi-squared = 781, df = 5, p-value <2e-166.3 Effect size
rstatix:::kruskal_effsize(salary ~ major,
data = df)# A tibble: 1 × 5
.y. n effsize method magnitude
* <chr> <int> <dbl> <chr> <ord>
1 salary 1753 0.444 eta2[H] large 6.4 Post-hoc analysis | Multiple comparisons
6.4.1 Bonferroni
• “bonferroni” Kiểm soát family-wise error rate (FWER), khá bảo thủ • “holm” Kiểm soát FWER, thường ít bảo thủ hơn Bonferroni • “bh” Kiểm soát false discovery rate (FDR)
library(FSA)
FSA:::dunnTest(salary ~ major,
data = df,
method = "bonferroni") Comparison Z P.unadj P.adj
1 Biology - Economics and Finance -2.14 3.2e-02 4.8e-01
2 Biology - Environmental Sciences 13.68 1.4e-42 2.1e-41
3 Economics and Finance - Environmental Sciences 15.45 7.7e-54 1.2e-52
4 Biology - Mathematics and Statistics -0.86 3.9e-01 1.0e+00
5 Economics and Finance - Mathematics and Statistics 1.24 2.1e-01 1.0e+00
6 Environmental Sciences - Mathematics and Statistics -14.06 6.6e-45 9.9e-44
7 Biology - Political Science 14.18 1.3e-45 1.9e-44
8 Economics and Finance - Political Science 15.88 9.3e-57 1.4e-55
9 Environmental Sciences - Political Science 1.12 2.6e-01 1.0e+00
10 Mathematics and Statistics - Political Science 14.55 6.1e-48 9.1e-47
11 Biology - Social Sciences 18.07 5.9e-73 8.8e-72
12 Economics and Finance - Social Sciences 19.56 3.4e-85 5.0e-84
13 Environmental Sciences - Social Sciences 6.26 3.9e-10 5.9e-09
14 Mathematics and Statistics - Social Sciences 18.34 4.2e-75 6.2e-74
15 Political Science - Social Sciences 5.07 4.0e-07 6.0e-066.5 Vẽ đồ thị
library(ggstatsplot)
library(tidyverse)
set.seed(20)
kq <- ggstatsplot:::ggbetweenstats(data = df,
x = major,
y = salary,
var.equal = TRUE,
type = "nonparametric",
plot.type = "box",
pairwise.comparisons = T,
pairwise.display = "all",
p.adjust.method = "bonferroni",
# p.adjust.method = "none",
centrality.plotting = T,
centrality.path = T,
digits = 4L,
centrality.point.args = list(size = 4, color = "darkred"),
centrality.label.args = list(size = 4,
nudge_x = 0.4,
segment.linetype = 4),
ggtheme = theme_classic(base_size = 12),
bf.message = T,
results.subtitle = T,
caption = "",
ggsignif.args = list(size = 0.5,
vjust = 0,
step_increase = 0.1,
textsize = 5),
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0, jitter.width = 0, jitter.height = 0),
alpha = 0.5,
size = 6,
stroke = 0,
na.rm = TRUE)
)
# kq
ggstatsplot:::extract_stats(kq)$subtitle_data
# A tibble: 1 × 15
parameter1 parameter2 statistic df.error p.value method effectsize estimate conf.level conf.low conf.high conf.method
<chr> <chr> <dbl> <int> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 salary major 781. 5 1.18e-166 Kruskal-Wallis rank sum test Epsilon2 (rank) 0.446 0.95 0.418 1 percentile bootstrap
conf.iterations n.obs expression
<int> <int> <list>
1 100 1753 <language>
$caption_data
NULL
$pairwise_comparisons_data
# A tibble: 15 × 9
group1 group2 statistic p.value alternative distribution p.adjust.method test expression
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <list>
1 Biology Economics and Finance 2.14 4.79e- 1 two.sided z Bonferroni Dunn <language>
2 Biology Environmental Sciences 13.7 2.06e-41 two.sided z Bonferroni Dunn <language>
3 Biology Mathematics and Statistics 0.857 1 e+ 0 two.sided z Bonferroni Dunn <language>
4 Biology Political Science 14.2 1.92e-44 two.sided z Bonferroni Dunn <language>
5 Biology Social Sciences 18.1 8.84e-72 two.sided z Bonferroni Dunn <language>
6 Economics and Finance Environmental Sciences 15.4 1.16e-52 two.sided z Bonferroni Dunn <language>
7 Economics and Finance Mathematics and Statistics 1.24 1 e+ 0 two.sided z Bonferroni Dunn <language>
8 Economics and Finance Political Science 15.9 1.39e-55 two.sided z Bonferroni Dunn <language>
9 Economics and Finance Social Sciences 19.6 5.03e-84 two.sided z Bonferroni Dunn <language>
10 Environmental Sciences Mathematics and Statistics 14.1 9.92e-44 two.sided z Bonferroni Dunn <language>
11 Environmental Sciences Political Science 1.12 1 e+ 0 two.sided z Bonferroni Dunn <language>
12 Environmental Sciences Social Sciences 6.26 5.85e- 9 two.sided z Bonferroni Dunn <language>
13 Mathematics and Statistics Political Science 14.5 9.12e-47 two.sided z Bonferroni Dunn <language>
14 Mathematics and Statistics Social Sciences 18.3 6.23e-74 two.sided z Bonferroni Dunn <language>
15 Political Science Social Sciences 5.07 6.05e- 6 two.sided z Bonferroni Dunn <language>
$descriptive_data
NULL
$one_sample_data
NULL
$tidy_data
NULL
$glance_data
NULL
attr(,"class")
[1] "ggstatsplot_stats" "list" ggsave(filename = "kq_dunn.png",
plot = kq,
width = 14,
height = 10,
units = "in",
dpi = 300)
6.6 Summary giữa các nhóm
library(DescTools)
DescTools:::MedianCI(upper_30)median lwr.ci upr.ci
43712 39356 51706
attr(,"conf.level")
[1] 0.97quantile(upper_30) 0% 25% 50% 75% 100%
22191 36805 43712 54950 71595 options(pillar.sigfig = 8)
df |> dplyr:::group_by(major) |>
dplyr:::summarise(median = DescTools:::MedianCI(salary)[1],
lower = DescTools:::MedianCI(salary)[2],
upper = DescTools:::MedianCI(salary)[3],
Q2 = quantile(salary)[2],
Q4 = quantile(salary)[4],
IQR = IQR(salary))# A tibble: 6 × 7
major median lower upper Q2 Q4 IQR
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Biology 47085.666 45676.770 48969.948 41201.484 53332.147 12130.664
2 Economics and Finance 49223.300 47842.436 50416.159 43166.767 55939.829 12773.062
3 Environmental Sciences 37691.928 36750.087 38476.916 33365.419 41924.629 8559.2101
4 Mathematics and Statistics 47670.778 46306.616 49317.398 42758.735 53590.981 10832.246
5 Political Science 36279.911 34444.640 37529.564 29934.009 42064.680 12130.670
6 Social Sciences 32031.546 31325.999 33200.275 27612.823 36021.198 8408.3754Các trường hợp phân tích paired samples ghép cặp before-after hay theo dõi theo thời gian được mình trình bày chi tiết ở link này [Paired samples]
