具有块设计和重复测量的 ANOVA

Question

我正在尝试对在同一生长季节在 2 个地点进行的田间试验进行一些统计分析。

在两个站点 (Site，级别：HF|NW)，实验设计是一个 RCBD，其中包含 4 个 (n=4) 块（Block，级别：1|2|3|4 在每个 Site 内）。 有 4 种处理——3 种不同形式的氮肥和一个对照（无氮肥）（Treatment，水平：AN、U、IU、C）。 在田间试验期间，有 3 个不同的时期，从添加肥料开始到收割草结束。这些时期在因子N_app 下被赋予了级别 1|2|3。

我想在一系列测量值上检验以下零假设 H0：

Treatment (H0) 对测量没有影响

我特别感兴趣的两个指标是：草产量和氨排放量。

从草产量 (Dry_tonnes_ha) 开始 显示here, a nice balanced data set

可以使用以下代码在 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")))

我已经尝试使用以下方法对此进行方差分析：

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 个不同时期从同一个地块进行的 - 将这种重复测量构建到模型中的最佳方法是什么？

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

我有another data set for NH3 emissions here。 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_app 1 n=3，但对于N_app 2 & 3 n=4 对于所有N_app 级别，NW n=4。 在NF，仅在Treatment 水平U 和IU 上进行测量 在NW 测量是在Treatment 水平AN、U 和IU 上进行的

我不确定如何处理这种增加的复杂性。我很想将其分析为 2 个单独的站点（每个站点的 N_app 句点不相同的事实可能会鼓励这种方法）。 我可以在这里使用类型 iii 的平方和方差分析吗？

有人向我建议，线性混合建模方法可能是前进的方向，但我不熟悉使用这些方法。

欢迎您对上述任何内容提出意见。感谢您的宝贵时间。

罗里

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 开始，就好像它们是独立的因素一样。但是，如果您知道它们不是独立的，则可能应该嵌套它们。

我认为您的问题和疑虑非常开放。我的建议是，只要您有合理的理由，就继续前进。

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 年）。我不确定在这种情况下是否需要它们。 Vignette provided

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)

在这些情节中我没有看到任何太麻烦的地方。