R：数字列的单向方差分析和成对 post 临时测试（Turkey、Scheffe 或其他）

Question

我在数据框中有三列 沙丘（页面底部下方） 描述了三种不同沙丘生态系统记录的马拉姆草覆盖率百分比：

(1) 已恢复； (2) 退化；和 (3) 天然的；

我进行了两种不同的单向方差分析测试（如下）——测试 1 和测试 2——以确定生态系统之间的显着差异。测试 1 清楚地显示了生态系统之间的显着差异；然而，检验 2 没有显着差异。箱形图（下图）显示了生态系统之间方差的明显差异。

之后，我融化了数据框以生成阶乘列（即，标题为 Ecosystem.Type），这也是响应变量。这个想法是应用一个 glm 模型（下面的测试 3）来测试单向方差分析；但是，此方法不成功（请在下面找到错误消息）。

问题

我很困惑我执行每个单向方差分析测试的代码是否正确以及执行 post 临时测试（土耳其 HSD、Scheffe 或其他）以区分显着不同的生态系统对的正确程序.如果有人可以提供帮助，我将非常感谢您的建议。非常感谢....

data(dune)

测试 1

dune.type.1<-aov(Natural~Restored+Degraded, data=dune)
summary.aov(dune.type.1, intercept=T)

               Df Sum Sq Mean Sq F value   Pr(>F)    
     (Intercept)  1  34694   34694 138.679 1.34e-09 ***
     Restored     1     94      94   0.375    0.548    
     Degraded     1    486     486   1.942    0.181    
     Residuals   17   4253     250                     
           ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Post-hoc 测试

    posthoc<-TukeyHSD(dune.type.1, conf.level=0.95)

    Error in TukeyHSD.aov(dune.type.1, conf.level = 0.95) : 

    no factors in the fitted model

    In addition: Warning messages:
    1: In replications(paste("~", xx), data = mf) :
       non-factors ignored: Restored
    2: In replications(paste("~", xx), data = mf) :
       non-factors ignored: Degraded

测试 2

     dune1<-aov(Restored~Natural, data=dune)
     dune2<-aov(Restored~Degraded, data=dune)
     dune3<-aov(Degraded~Natural, data=dune)

     summary(dune1)

                 Df Sum Sq Mean Sq F value Pr(>F)
     Natural      1     86   85.58   0.356  0.558
     Residuals   18   4325  240.26               

    summary(dune2)

                 Df Sum Sq Mean Sq F value Pr(>F)
     Degraded     1    160   159.7   0.676  0.422
     Residuals   18   4250   236.1               

     summary(dune3)

                 Df Sum Sq Mean Sq F value Pr(>F)
     Natural      1  168.5  168.49   2.318  0.145
     Residuals   18 1308.5   72.69

测试 3

melt.dune<-melt(dune, measure.vars=c("Degraded", "Restored", "Natural"))


colnames(melt.dune)=c("Ecosystem.Type", "Percentage.cover")
melt.dune$Percentage.cover<-as.numeric(melt.dune$Percentage.cover)

glm.dune<-glm(Ecosystem.Type~Percentage.cover, data=melt.dune)
summary(glm.dune)

Error

glm.dune<-glm(Ecosystem.Type~Percentage.cover, data=melt.dune)
Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  : 
NA/NaN/Inf in 'y'
In addition: Warning messages:
1: In Ops.factor(y, mu) : ‘-’ not meaningful for factors
2: In Ops.factor(eta, offset) : ‘-’ not meaningful for factors
3: In Ops.factor(y, mu) : ‘-’ not meaningful for factors

融化的数据框

structure(list(Ecosystem.Type = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Degraded", "Restored", 
"Natural"), class = "factor"), Percentage.cover = c(12, 17, 21, 
11, 22, 16, 7, 9, 14, 2, 3, 15, 23, 4, 19, 36, 26, 4, 15, 23, 
38, 46, 65, 35, 54, 29, 48, 13, 19, 33, 37, 55, 11, 53, 13, 24, 
28, 44, 42, 39, 18, 61, 31, 46, 51, 51, 41, 44, 55, 47, 73, 43, 
25, 42, 21, 13, 65, 30, 47, 29)), row.names = c(NA, -60L), .Names =         c("Ecosystem.Type", 
 "Percentage.cover"), class = "data.frame")

数据

 structure(list(Degraded = c(12L, 17L, 21L, 11L, 22L, 16L, 7L, 
 9L, 14L, 2L, 3L, 15L, 23L, 4L, 19L, 36L, 26L, 4L, 15L, 23L), 
 Restored = c(38L, 46L, 65L, 35L, 54L, 29L, 48L, 13L, 19L, 
 33L, 37L, 55L, 11L, 53L, 13L, 24L, 28L, 44L, 42L, 39L), Natural = c(18L, 
 61L, 31L, 46L, 51L, 51L, 41L, 44L, 55L, 47L, 73L, 43L, 25L, 
 42L, 21L, 13L, 65L, 30L, 47L, 29L)), .Names = c("Degraded", 
 "Restored", "Natural"), class = "data.frame", row.names = c(NA, 
 -20L))

Answer 1

我想向您指出几件事。

首先，测试 1 和测试 2 产生相似的结果。唯一的区别是您在测试 1 中选择了一个截距，因此结果告诉您，如果您拟合线性模型（我将在几分钟内谈到），则需要截距。因此，您看到的重要性在于您强制拟合的线是否需要截距。尝试对其他结果使用 "intercept=T"，我很确定您会得到类似的结果。

您应该注意的第二件事是您尝试拟合的线性模型。 dune.type.1 模型是您实际看到不同沙丘生态系统之间关联程度的模型。换句话说，您假设自然和恢复之间存在线性关联，并且随着恢复的每增加一个单位，自然也会有所增加。如果我理解正确的话，你想要的是检查平均差异而不是它们的相关性。因此你可以做两件事：

数据已准备好执行 t.tests（比较几个类别之间的平均值的测试）。这在 R 中很容易做到并且有效，因为所有变量都相当正常。然而，您将遇到多个测试问题（您可能会执行 3 次 t 测试以获得所有平均差异），因此需要使用 Bonferroni 校正。
但我认为你真正想要的是以下内容：

先改造数据

       data <- data.frame(v = c(dune$Degraded, dune$Restored, dune$Natural), 
                   labels = c(rep("Degraded", times = 20), rep("Restored", times = 20), 
                              rep("Natural", times = 20)))

然后拟合一个线性模型

    mod.1 <- lm(v ~ labels, data = data)
    summary(mod.1)
    lm(formula = v ~ labels, data = data)

Residuals:
Min      1Q  Median      3Q     Max 
-28.650 -10.725   0.875   8.050  31.350 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      14.950      3.066   4.875 9.07e-06 ***
labelsNatural    26.700      4.337   6.157 7.95e-08 ***
labelsRestored   21.350      4.337   4.923 7.64e-06 ***

你实际上可以看到基线类别（即退化）的平均值明显小于自然类别等的平均值。你还可以检查模型假设，看看你的模型是否是很合身

    qqnorm(residuals(mod.1))
    qqline(residuals(mod.1))

它们的残差相当正常，因此模型很好。您还可以按照方差分析方法进行操作并获得：

    anova.model <- aov(v ~ labels, data = data))
    summary(anova.model)

             Df Sum Sq Mean Sq F value   Pr(>F)    
 labels       2   7982    3991   21.22 1.29e-07 ***
 Residuals   57  10720     188

这表明沙丘生态系统的平均值之间至少存在一个显着差异，并针对逐点间隔跟进Tukey：

    post <- TukeyHSD(aov(v ~ labels, data = data))
    plot(post, ylim = c(0, 4))

已经针对多重测试进行了调整:)

R：数字列的单向方差分析和成对 post 临时测试（Turkey、Scheffe 或其他）

R: One Way Anova and pairwise post hoc tests (Turkey, Scheffe or other) for numerical columns

r

anova

posthoc