2

我有一组从 RCBD 设计中获得的数据。数据是莴苣中烫伤的发生率(即由环境因素引起的生理性叶片疾病,这是不好的)。

我的实验包括来自重组近交系(即 RIL)种群的 3 个块和 92 个莴苣品种。

我的目标是分析数据,找出与其他品种相比,哪些品种的烫伤发生率在统计学上有显着差异。

即使在转换后,我的数据也不符合方差同质性的 ANOVA 假设;因此,我的下一个方法是使用非参数测试。我发现我可以使用弗里德曼检验来分析数据。为了进行此测试并测试哪些品种与其他品种有统计差异,我在以下网站上方便地找到了一个功能:

http://www.r-statistics.com/2010/02/post-hoc-analysis-for-friedmans-test-r-code/

但是,按照说明格式化数据后,我无法运行该功能,如下所示:

格式化数据

    > tip.data.2011 = read.csv("Salinas_2011_tipburn_analysis.csv", header = TRUE)

    > head(tip.data.2011)
      RIL Block Tipburn_percentage
    1 110     1                0.0
    2 110     2                0.0
    3 110     3                0.0
    4 111     1               37.5
    5 111     2               12.5
    6 111     3               37.5

    > tip.data.2011.formated = data.frame( 
    +   Tipburn = tip.data.2011$Tipburn_percentage, 
    +     RIL = factor(rep(subset(tip.data.2011, Block == 1)[,1], rep(3, 92))),
    +   Block = factor(rep(1:3, 92))
    + )

    > head(tip.data.2011.formated)
      Tipburn RIL Block
    1     0.0 110     1
    2     0.0 110     2
    3     0.0 110     3
    4    37.5 111     1
    5    12.5 111     2
    6    37.5 111     3

运行函数

    > friedman.test.with.post.hoc(Tipburn ~ RIL | Block, tip.data.2011.formated)
    Error in mvt(lower = lower, upper = upper, df = 0, corr = corr, delta = mean,  : 
      only dimensions 1 <= n <= 1000 allowed

我收到错误消息“mvt 中的错误(下 = 下,上 = 上,df = 0,corr = corr,delta = mean,:只允许尺寸 1 <= n <= 1000”

这个错误信息是什么意思?

关于如何修复它的任何想法?

参考@DWin 的请求:

    > with(tip.data.2011.formated, tapply(Tipburn, list(RIL, Block), length))
        1 2 3
    110 1 1 1
    111 1 1 1
    112 1 1 1
    113 1 1 1
    114 1 1 1
    115 1 1 1
    116 1 1 1
    117 1 1 1
    118 1 1 1
    119 1 1 1
    120 1 1 1 ... etc.
    198 1 1 1
    199 1 1 1
    200 1 1 1
    SAL 1 1 1

    > str(tip.data.2011.formated)
    'data.frame':   276 obs. of  3 variables:
     $ Tipburn: num  0 0 0 37.5 12.5 37.5 0 0 12.5 75 ...
     $ RIL    : Factor w/ 92 levels "110","111","112",..: 1 1 1 2 2 2 3 3 3 4 ...
     $ Block  : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...

我对一个可以工作的数据集执行了与上面所示相同的过程(数据集在上面提到的网站中给出)。

    > WineTasting <- data.frame(
    +     Taste = c(5.40, 5.50, 5.55,
    +               5.85, 5.70, 5.75,
    +               5.20, 5.60, 5.50,
    +               5.55, 5.50, 5.40,
    +               5.90, 5.85, 5.70,
    +               5.45, 5.55, 5.60,
    +               5.40, 5.40, 5.35,
    +               5.45, 5.50, 5.35,
    +               5.25, 5.15, 5.00,
    +               5.85, 5.80, 5.70,
    +               5.25, 5.20, 5.10,
    +               5.65, 5.55, 5.45,
    +               5.60, 5.35, 5.45,
    +               5.05, 5.00, 4.95,
    +               5.50, 5.50, 5.40,
    +               5.45, 5.55, 5.50,
    +               5.55, 5.55, 5.35,
    +               5.45, 5.50, 5.55,
    +               5.50, 5.45, 5.25,
    +               5.65, 5.60, 5.40,
    +               5.70, 5.65, 5.55,
    +               6.30, 6.30, 6.25),
    +     Wine = factor(rep(c("Wine A", "Wine B", "Wine C"), 22)),
    +     Taster = factor(rep(1:22, rep(3, 22))))

    > with(WineTasting, tapply(Taste, list(Wine, Taster), length))
           1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
    Wine A 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1
    Wine B 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1
    Wine C 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1

    > str(WineTasting)
    'data.frame':   66 obs. of  3 variables:
     $ Taste : num  5.4 5.5 5.55 5.85 5.7 5.75 5.2 5.6 5.5 5.55 ...
     $ Wine  : Factor w/ 3 levels "Wine A","Wine B",..: 1 2 3 1 2 3 1 2 3 1 ...
     $ Taster: Factor w/ 22 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...

我的数据的格式似乎与运行的示例数据集相同。我不认为有大量的 1 观察到 Tipburn 是问题。

谢谢,

米格尔

4

1 回答 1

2

(我对 CrossValidated 交流中的评论不以为然。)问题是为什么您会收到 Tal Galili 的代码错误。从该网页下载代码表明错误消息必须来自它所依赖的包之一,因为该错误不会在该函数中引发。引发错误的函数是mvt,但该函数不在代码中提到的任何包中,至少当我查看这些包的最新版本时。(原来它是mvtnorm由 附加的multcomp。)mvt引发错误的代码是

if (n > 1000) 
    stop("only dimensions 1 <= n <= 1000 allowed")

所以我怀疑你的问题比作者预期的要multcompmvtnorm。你能看看:

  with(tip.data.2011.formated,  
               tapply(Tipburn, list(RIL,  Block), 
                                  length) ) )

鉴于您对该诊断查询的结果,我猜您的格式错误RILBlock变量。您有大量类别,其中Tipburn. 也尝试发布...这次作为对您的问题的编辑...的结果str(tip.data.2011.formated)

在使用少量类别“Wine”的因子成功的代码中,“主效应”变量是“主效应”变量。在您的设置中,将具有 92 个水平的因子作为“主效应”变量,将具有少量类别的变量作为分层变量。我认为您从未表达过正在测试的特定假设,特别是“块”值的含义是什么,所以我不确定问题是什么。如果颠倒这些变量的顺序是有意义的,即如果问题是块差异是否重要,您应该尝试以下公式:Tipburn ~ Block | RIL

我还担心您可能已将此数据集从原始数据减少到经过处理的版本,并且glm对原始计数数据的分析可能更有意义。分析比例表明这始于“计数数据”。

于 2013-06-27T04:07:44.810 回答