1

相关问题的链接（如何从数据子集随机抽取并在 R 中引导统计测试）提供了一个很好的示例，说明如何对数据帧中随机抽取的数据子样本进行统计测试。作为这个问题的扩展，我想知道如何对统计测试的引导迭代执行事后测试，在该测试中发现组之间存在显着差异。

假设我在三年内（Y1、Y2、Y3）对植物进行了采样。我想知道使用 Kruskal-Wallis 检验，植物的中位长度在不同年份之间是否存在显着差异。如果他们这样做（即 p 值 <0.05），我想知道哪些年份显示出显着差异，使用 Wilcoxon 秩和检验。由于我的数据框中有一些植物在某一年内进行了多次测量，因此我将在每一年内为这些植物随机绘制一行，用于每次统计测试迭代以防止伪复制。该过程将重复 10 次。

示例数据：

structure(list(Plant = c(1L, 2L, 3L, 4L, 5L, 6L, 6L, 7L, 8L, 
9L, 10L, 10L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 18L, 
19L, 20L, 21L), Year = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), Length = c(110L, 99L, 124L, 154L, 112L, 129L, 93L, 132L, 178L, 
206L, 177L, 257L, 173L, 222L, 167L, 192L, 354L, 299L, 265L, 301L, 
341L, 316L, 289L, 267L, 250L)), .Names = c("Plant", "Year", "Length"
), class = "data.frame", row.names = c(NA, -25L))

   Plant Year Length
1      1    1    110
2      2    1     99
3      3    1    124
4      4    1    154
5      5    1    112
6      6    1    129
7      6    1     93    etc….

我的问题是，在存在显着差异的情况下，如何在引导重复中执行事后测试（然后将每个测试统计量、p 值和参数值保存在矩阵/数据框中）。我已经尝试了下面的代码，但我得到的只是一个与迭代次数相同长度的输出矩阵，这不应该是这种情况。

library(plyr)

# function to draw sample for each iteration, perform the KW test and perform post-hoc tests when differences are significant
myrandomph <- function(P,Q){
  ss <- ddply(P, .(Plant), function(x) { y <- x[sample(nrow(x), 1) ,] })
  kw <- kruskal.test(ss$Length, ss$Year)

  if(kw$p.value < 0.05){ 
    w1 <- wilcox.test(ss$Length[ss$Year =="1"], ss$Length[ss$Year =="2"], paired=FALSE)
    return(c(stat = w1$statistic, p = w1$p.value))
    w2 <- wilcox.test(ss$Length[ss$Year =="1"], ss$Length[ss$Year =="3"], paired=FALSE)
    return(c(stat = w2$statistic, p = w2$p.value))
    w3 <- wilcox.test(ss$Length[ss$Year =="2"], ss$Length[ss$Year =="3"], paired=FALSE)
    return(c(stat = w3$statistic, p = w3$p.value))
  }else{  
  return(c(stat = kw$statistic, p = kw$p.value, df = kw$parameter)) }
}

# repeat myrandomph 10 times
test_results <- do.call( rbind, replicate(10, myrandomph(df), simplify=FALSE ) )
colnames(test_results) <- c("Statistic", "P.value", "Parameter")

在 Kruskal-Wallis 测试很重要的情况下，我希望得到一个带有 KW 测试输出的数据框行，以及事后测试的每个测试输出的一行（即具有统计值的列，一列使用 p 值和行标签指定运行了哪个事后测试：w1、w2 或 w3）。在 Kruskal-Wallis 检验不显着的情况下，我只希望返回 KW 统计数据、p 值和参数。任何建议将不胜感激！

4

0 回答 0