1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
| library(multcomp)
library(pgirmess)
stat_xiang <- function(df, value, group, method, level){
df_sub <- df[,c(value, group)]
colnames(df_sub) <- c('value','group')
number <- as.data.frame(table(df_sub$group))
colnames(number) <- c('group','number')
df_sub <- merge(df_sub, number, by = 'group')
mean_value <- aggregate(df_sub$value, by = list(df_sub$group), FUN = mean)
colnames(mean_value) <- c('group','mean')
sd_value <- aggregate(df_sub$value, by = list(df_sub$group), FUN = sd)
colnames(sd_value) <- c('group','standard_deviation')
temp_df <- merge(mean_value, sd_value, by = 'group')
df_sub <- merge(df_sub, temp_df, by = 'group')
df_sub$standard_error <- df_sub$standard_deviation / sqrt(df_sub$number)
if (length(unique(df_sub$group)) == 2) {
if (method == 't.test') {
fit <- t.test(value ~ group, data = df_sub)
pvalue <- fit[["p.value"]]
signif <- ifelse(pvalue < 0.001,'***',
ifelse(pvalue > 0.001 & pvalue < 0.01, '**',
ifelse(pvalue > 0.05, 'NS','*')))
}
if (method == 'wilcox') {
fit <- wilcox.test(value ~ group, data = df_sub)
pvalue <- fit[["p.value"]]
signif <- ifelse(pvalue < 0.001,'***',
ifelse(pvalue > 0.001 & pvalue < 0.01, '**',
ifelse(pvalue > 0.05, 'NS','*')))
}
sig <- data.frame(group = unique(df_sub$group),
method = method,
level = level,
pvalue = pvalue,
signif = c(signif,''))
}
if (length(unique(df_sub$group)) > 2) {
if (method == 'anova') {
fit <- aov(value ~ group, data = df_sub)
pvalue <- summary(fit)[[1]][["Pr(>F)"]][1]
tuk <- glht(fit, linfct = mcp(group = 'Tukey'))
signif <- cld(tuk, level = level, ddecreasing = TRUE)[["mcletters"]][["Letters"]]
signif <- as.data.frame(signif)
colnames(signif) = 'signif'
signif$group <- rownames(signif)
signif$method <- method
signif$pvalue <- pvalue
signif$level <- level
sig <- signif[,c('group','method','level','pvalue','signif')]
}
if (method == 'kruskal') {
fit <- kruskal.test(value ~ group, data = df_sub)
pvalue <- fit[["p.value"]]
if (pvalue < 0.05) {
fit_2 <- as.data.frame(kruskalmc(df_sub$value, df_sub$group, probs = 1-level))
signif <- as.data.frame(fit_2)
signif$statistic <- rownames(signif)
colnames(signif)[2] <- 'group_comp'
signif$group <- unique(df_sub$group)
sig <- data.frame(group = unique(df_sub$group),
method = method,
level = level,
pvalue = pvalue)
sig <- merge(signif, sig, by = 'group')
}else{
sig <- data.frame(group = unique(df_sub$group),
method = method,
level = level,
pvalue = pvalue,
signif = 'NS')
}
}
}
results <- merge(df_sub,sig, by = 'group', all.x = TRUE)
return(results)
}
|