博文
方差分析多重比较出图
||
下面这些代码将展示由Francois Gillet编写的两个方差分析多重比较的函数 boxplert()和boxplerk()，来源Numerical Ecology with R (second Edition)
#普通方差分析多重比较
boxplert <- function(X,
Y,
main = NULL,
xlab = NULL,
ylab = NULL,
bcol = "bisque",
p.adj = "none",
cexy = 1,
varwidth = TRUE,
las = 1,
paired = FALSE)
{#Y=factor(Y,levels=unique(Y)) #如果运行次代码，表示组别是按照原来的顺序，而不是字母顺序排列
aa <- levels(as.factor(Y))
an <- as.character(c(1:length(aa)))
tt1 <- matrix(nrow = length(aa), ncol = 6)
for (i in 1:length(aa))
{
temp <- X[Y == aa[i]]
tt1[i, 1] <- mean(temp, na.rm = TRUE)
tt1[i, 2] <- sd(temp, na.rm = TRUE) / sqrt(length(temp))
tt1[i, 3] <- sd(temp, na.rm = TRUE)
tt1[i, 4] <- min(temp, na.rm = TRUE)
tt1[i, 5] <- max(temp, na.rm = TRUE)
tt1[i, 6] <- length(temp)
}
tt1 <- as.data.frame(tt1)
row.names(tt1) <- aa
colnames(tt1) <- c("mean", "se", "sd", "min", "max", "n")
boxplot(
X ~ Y,
main = main,
xlab = xlab,
ylab = ylab,
las = las,
col = bcol,
cex.axis = cexy,
cex.lab = cexy,
varwidth = varwidth
)
require(agricolae)
Yn <- factor(Y, labels = an)
sig <- "ns"
model <- aov(X ~ Yn)
if (paired == TRUE & length(aa) == 2)
{
coms <- t.test(X ~ Yn, paired = TRUE)
pp <- coms$p.value
} else
{
pp <- anova(model)$Pr[1]
}
if (pp <= 0.1)
sig <- "."
if (pp <= 0.05)
sig <- "*"
if (pp <= 0.01)
sig <- "**"
if (pp <= 0.001)
sig <- "***"
mtext(
sig,
side = 3,
line = 0.5,
adj = 0,
cex = 2,
font = 1
)
if (pp <= 0.05) {
comp <- LSD.test(model, "Yn",
alpha = 0.05,
p.adj = p.adj,
group = TRUE) # gror <- comp$groups[order(comp$groups$groups), ]
# tt1$cld <- gror$M
gror <- comp$groups[order(rownames(comp$groups)), ]
tt1$group <- gror$groups
mtext(
tt1$group,
side = 3,
at = c(1:length(aa)),
line = 0.5,
cex = 1,
font = 4
)
}
list(comparison = tt1, p.value = pp)
}##非参数组间差异的Kruskal-Wallis检验
boxplerk <- function(X,
Y,
main = NULL,
xlab = NULL,
ylab = NULL,
bcol = "bisque",
p.adj = "none",
cexy = 1,
varwidth = TRUE,
las = 1,
paired = FALSE)
{#Y=factor(Y,levels=unique(Y)) #如果运行次代码，表示组别是按照原来的顺序，而不是字母顺序排列
aa <- levels(as.factor(Y))
an <- as.character(c(1:length(aa)))
tt1 <- matrix(nrow = length(aa), ncol = 7)
for (i in 1:length(aa))
{
temp <- X[Y == aa[i]]
tt1[i, 1] <- mean(temp, na.rm = TRUE)
tt1[i, 2] <- sd(temp, na.rm = TRUE) / sqrt(length(temp))
tt1[i, 3] <- sd(temp, na.rm = TRUE)
tt1[i, 4] <- min(temp, na.rm = TRUE)
tt1[i, 5] <- max(temp, na.rm = TRUE)
tt1[i, 6] <- median(temp, na.rm = TRUE)
tt1[i, 7] <- length(temp)
}
tt1 <- as.data.frame(tt1)
row.names(tt1) <- aa
colnames(tt1) <- c("mean", "se", "sd", "min", "max", "median", "n")
boxplot(
X ~ Y,
main = main,
xlab = xlab,
ylab = ylab,
las = las,
col = bcol,
cex.axis = cexy,
cex.lab = cexy,
varwidth = varwidth
)
require(agricolae)
Yn <- factor(Y, labels = an)
comp <- kruskal(X, Yn, p.adj = p.adj)
sig <- "ns"
if (paired == TRUE & length(aa) == 2)
{
coms <- wilcox.test(X ~ Yn, paired = TRUE)
pp <- coms$p.value
} else
{
pp <- comp$statistics$p.chisq
}
if (pp <= 0.1)
sig <- "."
if (pp <= 0.05)
sig <- "*"
if (pp <= 0.01)
sig <- "**"
if (pp <= 0.001)
sig <- "***"
gror <- comp$groups[order(rownames(comp$groups)), ]
tt1$rank <- gror$X
tt1$group <- gror$groups
mtext(
sig,
side = 3,
line = 0.5,
adj = 0,
cex = 2,
font = 1
)
if (pp <= 0.1)
mtext(
tt1$group,
side = 3,
at = c(1:length(aa)),
line = 0.5,
cex = 1,
font = 4
)
list(comparison = tt1, p.value = pp)
}
# # 普通方差分析的多重比较
library(agricolae)## Warning: package 'agricolae' was built under R version 3.4.4
data(sweetpotato)# # 检验方差分析假设 shapiro.test(resid(aov(sweetpotato$yield ~ sweetpotato$virus)))
## ## Shapiro-Wilk normality test ## ## data: resid(aov(sweetpotato$yield ~ sweetpotato$virus)) ## W = 0.95878, p-value = 0.7663
bartlett.test(sweetpotato$yield, sweetpotato$virus)
## ## Bartlett test of homogeneity of variances ## ## data: sweetpotato$yield and sweetpotato$virus ## Bartlett's K-squared = 2.3886, df = 3, p-value = 0.4958
boxplert( sweetpotato$yield, sweetpotato$virus, ylab = "yield", xlab = "virus", bcol = "orange", p.adj = "holm" )
## $comparison ## mean se sd min max n group ## cc 24.40000 2.084067 3.609709 21.7 28.5 3 b ## fc 12.86667 1.246774 2.159475 10.6 14.9 3 c ## ff 36.33333 4.233727 7.333030 28.0 41.8 3 a ## oo 36.90000 2.482606 4.300000 32.1 40.4 3 a ## ## $p.value ## [1] 0.0007334172
非参数Kruskal-Wallis检验多重比较
# # 案例 library(agricolae) library(stats) data(InsectSprays) boxplerk( InsectSprays$count, InsectSprays$spray, ylab = "count", xlab = "spray", bcol = "bisque", p.adj = "holm", las = 1)
## $comparison ## mean se sd min max median n rank group ## A 14.500000 1.3623732 4.719399 7 23 14.0 12 52.16667 a ## B 15.333333 1.2329647 4.271115 7 21 16.5 12 54.83333 a ## C 2.083333 0.5701984 1.975225 0 7 1.5 12 11.45833 c ## D 4.916667 0.7225621 2.503028 2 12 5.0 12 25.58333 b ## E 3.500000 0.5000000 1.732051 1 6 3.0 12 19.33333 bc ## F 16.666667 1.7936476 6.213378 9 26 15.0 12 55.62500 a ## ## $p.value ## [1] 1.510845e-10
请大家这里下载R markdown的文件
https://blog.sciencenet.cn/blog-267448-1149684.html
上一篇：NMDS的是否有解释率？
下一篇：R语言在生态学研究当中的应用分析
