带有区组设计和重复测量的ANOVA

Question

我试图对同一生长季内在两个场地进行的田间试验进行一些统计分析。

在两个场地（Site，水平：HF | NW）上，实验设计为RCBD，每个Site内有4个块（Block，水平：1 | 2 | 3 | 4）。共有4种处理方法-3种不同形式的氮肥和一种对照（无氮肥）（Treatment，水平：AN，U，IU，C）。在田间试验期间，有3个明显的时期，始于施肥，以草的收获结束。这些时期已在因子N_app下被赋予水平1 | 2 | 3。

我想测试以下零假设H0的一系列测量结果：

处理（Treatment）（H0）对测量结果没有影响

我特别感兴趣的是两项测量结果：草产量和氨排放。

从这里一个好的平衡数据集开始，以草产量（Dry_tonnes_ha）为例。

可以使用以下代码在R中下载数据：

library(tidyverse)

download.file('https://www.dropbox.com/s/w5ramntwdgpn0e3/HF_NW_grass_yield_data.csv?raw=1', destfile = "HF_NW_grass_yield_data.csv", method = "auto")
raw_data <- read.csv("HF_NW_grass_yield_data.csv", stringsAsFactors = FALSE)

HF_NW_grass <- raw_data %>% mutate_at(vars(Site, N_app, Block, Plot, Treatment), as.factor) %>% 
  mutate(Date = as.Date(Date, format = "%d/%m/%Y"),
         Treatment = factor(Treatment, levels = c("AN", "U", "IU", "C")))

我尝试使用以下方法运行ANOVA来分析这个问题：

model_1 <- aov(formula = Dry_tonnes_ha ~ Treatment * N_app + Site/Block, data = HF_NW_grass, projections = TRUE)

我对此有一些疑虑。

首先，测试假设的最佳方法是什么？对于简单的单向方差分析，我会在因变量（Dry_tonnes_ha）上使用shapiro.test()和bartlett.test()来评估正态性和方差异质性。这里可以使用相同的方法吗？

其次，我担心N_app是重复测量，因为同一地块在3个不同的时期进行了相同的测量 - 最好的方法是将这种重复测量建立到模型中的什么位置？

第三，我不确定在Site中嵌套Block的最佳方法。在两个站点上，Block的级别都是1:4。我是否需要为每个站点设置唯一的Block级别？

我在这里另一个NH3排放数据集。下载R代码：

download.file('https://www.dropbox.com/s/0ax16x95m2z3fb5/HF_NW_NH3_emissions.csv?raw=1', destfile = "HF_NW_NH3_emissions.csv", method = "auto")
raw_data_1 <- read.csv("HF_NW_NH3_emissions.csv", stringsAsFactors = FALSE)

HF_NW_NH3 <- raw_data_1 %>% mutate_at(vars(Site, N_app, Block, Plot, Treatment), as.factor) %>% 
  mutate(Treatment = factor(Treatment, levels = c("AN", "U", "IU", "C")))

对于这个问题，我有所有上述的担忧，并且数据集不平衡。 在HF中，对于N_app1 n=3，但对于N_app2和3 n=4。 在NW中，所有N_app级别的n=4。 在NF中，仅对Treatment级别U和IU进行了测量。 在NW中，对Treatment级别AN，U和IU进行了测量。

我不确定如何处理这个额外的复杂性。我倾向于将其分析为2个单独的站点（每个站点的N_app周期不同可能会鼓励采用这种方法）。 我可以在这里使用类型III平方和方差分析吗？

有人建议我采用线性混合模型方法，但我不熟悉使用这些方法。

我期待您对以上任何问题的看法。谢谢您的时间。

Rory

- Rory Shaw

Answer 1

回答你关于测试假设的最佳方法的第一个问题。虽然你尝试使用R中实施的另一种统计检验是合理的，但我实际上会直接可视化分布并查看数据是否符合ANOVA假设。这种方法可能有点主观，但在大多数情况下确实有效。

独立同分布（i.i.d）数据：这是一个问题，你可能已经基于你对数据的了解有了答案。可以使用卡方检验来确定独立性（或否）。

正态分布数据：使用直方图/QQ图进行检查。基于分布，尽管略微呈双峰分布，但我认为使用aov是合理的。

(似乎对数转换有助于进一步满足正态性假设。这是你可以考虑的事情，特别是对于下游分析。)

par(mfrow=c(2,2))
plot(density(HF_NW_grass$Dry_tonnes_ha), col="red", main="Density")
qqnorm(HF_NW_grass$Dry_tonnes_ha, col="red", main="qqplot")
qqline(HF_NW_grass$Dry_tonnes_ha)

DTH_trans <- log10(HF_NW_grass$Dry_tonnes_ha)
plot(density(DTH_trans), col="blue", main="transformed density")
qqnorm(DTH_trans, col="blue", main="transformed density")
qqline(DTH_trans)

关于您的第二个问题，关于如何将重复测量建立到模型中的最佳方式是什么：不幸的是，很难确定这样的“最佳”模型，但基于我的知识（主要通过基因组大数据），您可能需要使用线性混合效应模型。例如，可以通过lme4 R软件包来实现。由于您似乎已经知道如何在R中构建线性模型，因此应该没有问题应用lme4函数。

关于您的第三个问题，关于是否嵌套两个变量有些棘手。如果我是您，我会将Site和Block视为独立因子开始。但是，如果您知道它们不是独立的，那么您可能应该将它们嵌套。

我认为您的问题和顾虑相当开放。我的建议是，只要您有一个合理的理由，就可以继续进行。

- David C.

---

感谢您的回复。统计数据越复杂，过程似乎就越主观！关于假设的信息很有用。我现在开始质疑是否需要重复测量。我分析的所有测量仅对每个“N_app”进行一次，相同的“plot”用于3个“N_app”期间，这让我认为需要进行重复测量。我将再次研究线性混合模型。我还可以简化并单独分析两个站点。 - Rory Shaw

我猜线性混合模型也可以帮助处理不平衡的数据集。虽然我猜在这里可以使用type iii ss的Anova? - Rory Shaw

Answer 2

我同意 @David C 使用可视化诊断。简单的 QQ 图应该可以起到作用。

# dependent variable.
par(mfrow=c(1,2))
qqnorm(dt[,dry_tonnes_ha]); qqline(dt[,dry_tonnes_ha], probs= c(0.15, 0.85))
qqnorm(log(dt[,dry_tonnes_ha])); qqline(log(dt[,dry_tonnes_ha]), probs= c(0.15, 0.85))

对我来说，对数变换看起来很合理。你也可以从密度图中看到这一点，它的尾部很长，有些双峰。

par(mfrow=c(1,1))
plot(density(dt[,dry_tonnes_ha]))

你可以选择使用阵型图（Buja等人，2009），如果需要的话。我不确定在这种情况下是否需要。提供的文献

library(nullabor)
# this may not be the best X variable. I'm not familiar with your data
dt_l <- lineup(null_permute("dry_tonnes_ha"), dt)
qplot(dry_tonnes_ha, treatment, data = dt_l) + facet_wrap(~ .sample)

对于其他假设，您可以使用lm中的标准诊断图。

lm2 <- lm(log(dry_tonnes_ha) ~ treatment * n_app + site/block, data = dt)
plot(lm2)

我认为这些图表并没有什么麻烦的地方。

- alexwhitworth

---

谢谢@Alex，我以前没见过这种线路图-拥有更多选项总是很有用的。 - Rory Shaw