# 裂区试验设计
## 加载环境与数据准备
### 安装包检测并自动安装（必须联网）
```{r}
pkgs <- c("tidyverse","agricolae","here","knitr","latex2exp","cowplot")
not_installed <- pkgs[!(pkgs %in% installed.packages()[ , "Package"])]
if(length(not_installed)) install.packages(not_installed)
```
### 加载R必备包
```{r}
#| warning: false
#| label: loadingPackage
library(tidyverse)
library(agricolae)
library(here)
library(knitr)
library(latex2exp)
library(cowplot)
options(digits = 4)
```
### 试验数据读取
```{r}
#| label: tbl-dataset
#| tbl-cap: "试验数据集"
dt <- readxl::read_excel(here("RES/dataset/产量干物质.xlsx"),
sheet = "干物质",
skip = 2)
dt <- dt %>% filter(Sampling == 1) %>% select(-Sampling)
dt
```
## 试验数据检测
详细过程见[@sec-datatest]部分。
## 方差分析
由此可见，试验数据满足正态和方差齐次性，可正常进行方差分析。
以[@tbl-dataset]数据集为例，采用随机区组试验设计模型对试验数据进行方差分析，结果如[@tbl-ANOVA]。
```{r}
#| label: tbl-ANOVA
#| tbl-cap: "随机区组试验方差分析表"
response <- dt$Grainyield
ANOVA <- aov(response ~ Rep + Irrigation * Biochar + Error(Rep/Irrigation), data = dt) %>%
summary()
ANOVA <- pmap_df(list(data=ANOVA,name=names(ANOVA)),~.x[[1]],.id = "error") %>%
mutate(error=error %>% as.character())
ANOVA %>% select(-error) %>%
kable(digits=3)
```
有的学者可能对方差贡献率感兴趣，利用各因子及交互效应等平方和数据，可进一步得到各因子和交互效应贡献率情况，如[@tbl-ANOVA-Contribution]。
```{r}
#| label: tbl-ANOVA-Contribution
#| tbl-cap: "方差分析+方差贡献率"
ANOVA %>% select(-error) %>%
mutate(Contri=`Sum Sq`/sum(`Sum Sq`)*100) %>%
kable(digits = 3)
```
## 事后均值对比
### 定义HSD结果提取函数
```{r}
getHSD <- function(hsd){
merge(x=hsd$means %>% select(std),
y=hsd$groups,
by="row.names") %>%
rename(trt=Row.names)
}
```
### 对第1因子进行HSD均值检测
```{r}
#| label: tbl-hsd-irrigation
#| tbl-cap: "不同灌溉模式下产量均值对比结果"
# Irrigation 主区误差项
DFerror1 <- ANOVA %>% filter(error=="Error: Rep:Irrigation" & is.na(`F value`)) %>% getElement("Df")
MSerror1 <- ANOVA %>% filter(error=="Error: Rep:Irrigation" & is.na(`F value`)) %>% getElement("Mean Sq")
# Biochar和交互效应共一个误差项
DFerror2 <- ANOVA %>% filter(error=="Error: Within" & is.na(`F value`)) %>% getElement("Df")
MSerror2 <- ANOVA %>% filter(error=="Error: Within" & is.na(`F value`)) %>% getElement("Mean Sq")
means <- dt$Grainyield
hsd.Irrigation <- HSD.test(
means,
trt = dt$Irrigation,
DFerror = DFerror1,
MSerror = MSerror1
) %>%
getHSD() %>% select(trt, means, std, groups)
hsd.Irrigation %>% kable(digits=3)
```
### 对第2因子进行HSD均值检测
```{r}
#| label: tbl-hsd-biochar
#| tbl-cap: "不同生物炭水平下产量均值对比结果"
#|
hsd.Biochar <- HSD.test(means,
trt = dt$Biochar,
DFerror = DFerror2,
MSerror = MSerror2) %>%
getHSD() %>% select(trt, means, std, groups)
hsd.Biochar %>% kable(digits=3)
```
### 对交互因子进行HSD均值检测
```{r}
#| label: tbl-hsd-interaction
#| tbl-cap: "不同灌溉模式和生物炭水平下产量均值对比结果"
#|
hsd.Interaction <- HSD.test(
means,
trt = interaction(dt$Irrigation, dt$Biochar),
DFerror = DFerror2,
MSerror = MSerror2
) %>%
getHSD() %>% select(trt, means, std, groups) %>% separate(col = trt, into =
c("Irrigation", "Biochar"))
hsd.Interaction%>% kable(digits=3)
```
## 出图
### 设置试验因子各水平排序
```{r}
hsd.Irrigation$trt <- hsd.Irrigation$trt %>% fct_relevel("CF", "AWD")
hsd.Biochar$trt <- hsd.Biochar$trt %>% fct_relevel("B0", "B20", "B20M")
hsd.Interaction$Irrigation <- hsd.Interaction$Irrigation %>% fct_relevel("CF", "AWD")
hsd.Interaction$Biochar <- hsd.Interaction$Biochar %>% fct_relevel("B0", "B20", "B20M")
mylabels.irr <- c(CF = TeX("$I_{CF}$"), AWD = TeX("$I_{AWD}$"))
mylabels.bio <- c(
B0 = TeX("$B_{0}$"),
B20 = TeX("$B_{20}$"),
B20M = TeX("$B_{20M}$")
)
```
### 第1因子灌溉模式不同水平均值对比
```{r}
#| label: fig-ggplot-irrigation
#| fig-cap: "不同灌溉模式下产量情况"
G.irr <- ggplot(hsd.Irrigation ,aes(x=trt,y=means,fill=trt ))+
geom_bar(stat="identity")+
geom_errorbar(aes(ymin=means-std,ymax=means+std),
width=0.2)+
geom_text(aes(y=means+std,label=groups),vjust=-0.2)+
cowplot::theme_cowplot(font_size = 8,line_size = 0.4)+
scale_fill_discrete(labels=mylabels.irr )+
scale_x_discrete(labels=mylabels.irr )+
ylim(0,12)+
labs(x = "Irrigation Regime",
y = TeX("Grain yield ($t~ha^{-1}$)"),
fill = "Irrigation")
G.irr
```
### 第2因子生物炭不同水平均值对比
```{r}
#| label: fig-ggplot-biochar
#| fig-cap: "不同生物炭水平下产量情况"
G.Bio <- ggplot(hsd.Biochar ,aes(x=trt,y=means,fill=trt ))+
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = means - std, ymax = means + std), width = 0.2) +
geom_text(aes(y = means + std, label = groups), vjust = -0.2) +
cowplot::theme_cowplot(font_size = 8, line_size = 0.4) +
scale_fill_discrete(labels = mylabels.bio) +
scale_x_discrete(labels = mylabels.bio) +
ylim(0, 12) +
labs(
x = TeX("Biochar application rate ($t~ha^{-1}$)"),
y = TeX("Grain yield ($t~ha^{-1}$)"),
fill = "Biochar"
)
G.Bio
```
### 交互因子不同处理间平均值对比
```{r}
#| label: fig-ggplot-interaction
#| fig-cap: "不同灌溉模式和生物炭水平下产量情况"
dod <- position_dodge(width = 0.95)
G.int <- ggplot(hsd.Interaction , aes(x = Irrigation, y = means, fill = Biochar)) +
geom_bar(stat = "identity", position = dod) +
geom_errorbar(aes(ymin = means - std, ymax = means + std),
position = dod,
width = 0.2) +
geom_text(aes(y = means + std, label = groups),
vjust = -0.2,
position = dod) +
cowplot::theme_cowplot(font_size = 8, line_size = 0.4) +
scale_x_discrete(labels = mylabels.irr) +
scale_fill_discrete(labels = mylabels.bio) +
ylim(0, 12) +
labs(
x = TeX("Biochar application rate ($t~ha^{-1}$)"),
y = TeX("Grain yield ($t~ha^{-1}$)")
)
G.int
```
### 合成并出图
```{r}
#| label: fig-ggplot-merge
#| fig-cap: "主效应和交互效应"
GG <- plot_grid(
plot_grid(G.irr, G.Bio, nrow = 1, labels = letters),
G.int,
nrow = 2,
labels = c("", "c")
)
GG
```
