解析に必要なパッケージ
require(tidyverse)
require(plotly)解析事例
収量の分散分析（2024）
2024年度の水稲圃場実習では，要因として品種（4 品種）と施肥量（低リンと高リン）を取り上げた。なお、同様の実習が2019年度より行われているが、施肥は低窒素と高窒素の比較が行われていた。また、2022年以前はモミロマンを除いて異なる品種が栽培されていた。また、2022年度以降は水田の改修があったため影響で土壌の性質が変化している可能性もある。
施肥管理を効率的に行うため施肥管理を主区（1 次因子）、品種を副区（2 次因子）とした分割区法 （split-plot design）を用い、3 反復の実験を行った。
構造モデル
2023年度における試験区のある形質の測定値を \(Y_{ijk}\) とすると， \(Y_{ijk}\) の統計（構造）モデルは，1次因子水準数 \(a = 2\)，2次因子水準数 \(b = 4\)，反復 \(r = 3\)，として，
\[ Y_{ijk} = \mu + \alpha_i + \rho_k + e_{ik}^1 + \beta_j + (\alpha \beta)_{ij} + e_{ijk}^2, \ i=1,\ldots,a, j=1,\ldots,b, k=1,\ldots,r \] とおける．ただし，\(\mu\) は総平均，\(\alpha_i\) は施肥効果，\(\beta_j\) は品種効果，\((\alpha\beta)_{ij}\) は 施肥量と品種の交互作用，\(\rho_k\) はブロック効果，\(e_{ik}^1\) は 1次誤差、\(e_{ijk}^2\) は 2次誤差である．ブロック効果は、\(\rho_k \sim N(0, \sigma_B^2)\) を満たす変量（確率変数）で、誤差は，\(e_{ik}^1 \sim N(0, \sigma^2_1)\)、\(e^2_{ijk} \sim N(0, \sigma^2_2)\) を満たす互いに独立な変量であるとする。
1次誤差を用いて施肥効果やブロック効果を検定し、2次誤差を用いて品種効果や施肥と品種の交互作用を検定する。このため乱隗法に比べ、1次因子（施肥量）の効果に対する精度は落ちるが、2次因子（品種）や交互作用に対する精度は向上すると考えられる。
データの読み込みとグラフ描画 (2023)
データの読み込み
df <- read.csv("riceexp23.csv", stringsAsFactors = TRUE)
head(df)## Plot Year Rep Fertilizer Variety GrowSta_N_stems GrowSta_Plant_height
## 1 1 2023 1 Low_N Tachiayaka 251.8519 48.93333
## 2 2 2023 1 Low_N Momiroman 185.1852 31.43333
## 3 3 2023 1 Low_N Hokuriku193 266.6667 34.93333
## 4 4 2023 1 Low_N Koshihikari 207.4074 36.93333
## 5 5 2023 1 High_N Hokuriku193 281.4815 35.53333
## 6 6 2023 1 High_N Momiroman 392.5926 33.46667
## GrowSta_SPAD N_panicles N_grains_per_panicle N_grains_total Ripening_rate
## 1 43.96667 240.7407 164.2473 39541.03 0.6446330
## 2 35.53333 255.5556 258.2713 66002.67 0.3463521
## 3 38.50000 311.1111 136.1658 42362.69 0.8569052
## 4 38.23333 414.8148 126.1533 52330.27 0.5750007
## 5 38.46667 314.8148 169.6984 53423.57 0.6945010
## 6 37.33333 248.1481 257.3938 63871.80 0.2982613
## Thousand_grain_weight Grain_yield_gross Grain_yield Sink_capacity
## 1 28.00000 812.8922 713.7046 1107.149
## 2 27.05882 963.7093 544.6303 1785.955
## 3 27.29412 1053.5568 990.7985 1156.252
## 4 24.94118 933.6343 750.4786 1305.179
## 5 26.82353 1121.9526 795.3247 1433.009
## 6 27.52941 883.0191 524.4486 1758.353
## Sink_filling_rate Dry_weight_stems_leaves Harvest_index
## 1 0.7342213 884.4444 0.5195295
## 2 0.5396045 963.3333 0.5406374
## 3 0.9111825 1206.6667 0.5067069
## 4 0.7153307 894.4444 0.5511704
## 5 0.7829349 931.1111 0.5863669
## 6 0.5021853 971.1111 0.5168498データのサマリをみる。
summary(df)## Plot Year Rep Fertilizer Variety
## Min. : 1.00 Min. :2023 Min. :1 High_N:12 Hokuriku193:6
## 1st Qu.: 6.75 1st Qu.:2023 1st Qu.:1 Low_N :12 Koshihikari:6
## Median :12.50 Median :2023 Median :2 Momiroman :6
## Mean :12.50 Mean :2023 Mean :2 Tachiayaka :6
## 3rd Qu.:18.25 3rd Qu.:2023 3rd Qu.:3
## Max. :24.00 Max. :2023 Max. :3
## GrowSta_N_stems GrowSta_Plant_height GrowSta_SPAD N_panicles
## Min. :185.2 Min. :30.07 Min. :35.53 Min. :218.5
## 1st Qu.:246.3 1st Qu.:35.48 1st Qu.:38.41 1st Qu.:261.1
## Median :311.1 Median :41.38 Median :40.42 Median :277.8
## Mean :323.8 Mean :41.79 Mean :41.07 Mean :295.8
## 3rd Qu.:392.6 3rd Qu.:46.02 3rd Qu.:43.14 3rd Qu.:326.9
## Max. :555.6 Max. :57.60 Max. :51.60 Max. :414.8
## N_grains_per_panicle N_grains_total Ripening_rate Thousand_grain_weight
## Min. : 94.71 Min. :25257 Min. :0.2983 Min. :24.47
## 1st Qu.:134.30 1st Qu.:39101 1st Qu.:0.4334 1st Qu.:26.06
## Median :160.80 Median :49402 Median :0.6133 Median :26.82
## Mean :173.58 Mean :50028 Mean :0.5937 Mean :27.00
## 3rd Qu.:195.19 3rd Qu.:55651 3rd Qu.:0.6862 3rd Qu.:27.59
## Max. :317.76 Max. :88265 Max. :0.8637 Max. :32.47
## Grain_yield_gross Grain_yield Sink_capacity Sink_filling_rate
## Min. : 583.3 Min. :299.2 Min. : 689.4 Min. :0.4914
## 1st Qu.: 809.9 1st Qu.:523.4 1st Qu.:1134.9 1st Qu.:0.5867
## Median : 953.5 Median :716.7 Median :1282.7 Median :0.7408
## Mean : 944.3 Mean :653.5 Mean :1350.0 Mean :0.7274
## 3rd Qu.:1122.0 3rd Qu.:795.7 3rd Qu.:1588.8 3rd Qu.:0.8444
## Max. :1207.8 Max. :990.8 Max. :2409.1 Max. :0.9386
## Dry_weight_stems_leaves Harvest_index
## Min. : 713.3 Min. :0.2569
## 1st Qu.: 874.7 1st Qu.:0.5143
## Median : 951.1 Median :0.5444
## Mean : 996.0 Mean :0.5284
## 3rd Qu.:1043.3 3rd Qu.:0.5847
## Max. :1984.4 Max. :0.6348いくつかの変数が誤って数値データとして扱われている。それらを因子として変換を行う。
df <- df %>%
mutate(Plot = factor(Plot), Year = factor(Year), Rep = factor(Rep))summary(df)## Plot Year Rep Fertilizer Variety GrowSta_N_stems
## 1 : 1 2023:24 1:8 High_N:12 Hokuriku193:6 Min. :185.2
## 2 : 1 2:8 Low_N :12 Koshihikari:6 1st Qu.:246.3
## 3 : 1 3:8 Momiroman :6 Median :311.1
## 4 : 1 Tachiayaka :6 Mean :323.8
## 5 : 1 3rd Qu.:392.6
## 6 : 1 Max. :555.6
## (Other):18
## GrowSta_Plant_height GrowSta_SPAD N_panicles N_grains_per_panicle
## Min. :30.07 Min. :35.53 Min. :218.5 Min. : 94.71
## 1st Qu.:35.48 1st Qu.:38.41 1st Qu.:261.1 1st Qu.:134.30
## Median :41.38 Median :40.42 Median :277.8 Median :160.80
## Mean :41.79 Mean :41.07 Mean :295.8 Mean :173.58
## 3rd Qu.:46.02 3rd Qu.:43.14 3rd Qu.:326.9 3rd Qu.:195.19
## Max. :57.60 Max. :51.60 Max. :414.8 Max. :317.76
##
## N_grains_total Ripening_rate Thousand_grain_weight Grain_yield_gross
## Min. :25257 Min. :0.2983 Min. :24.47 Min. : 583.3
## 1st Qu.:39101 1st Qu.:0.4334 1st Qu.:26.06 1st Qu.: 809.9
## Median :49402 Median :0.6133 Median :26.82 Median : 953.5
## Mean :50028 Mean :0.5937 Mean :27.00 Mean : 944.3
## 3rd Qu.:55651 3rd Qu.:0.6862 3rd Qu.:27.59 3rd Qu.:1122.0
## Max. :88265 Max. :0.8637 Max. :32.47 Max. :1207.8
##
## Grain_yield Sink_capacity Sink_filling_rate Dry_weight_stems_leaves
## Min. :299.2 Min. : 689.4 Min. :0.4914 Min. : 713.3
## 1st Qu.:523.4 1st Qu.:1134.9 1st Qu.:0.5867 1st Qu.: 874.7
## Median :716.7 Median :1282.7 Median :0.7408 Median : 951.1
## Mean :653.5 Mean :1350.0 Mean :0.7274 Mean : 996.0
## 3rd Qu.:795.7 3rd Qu.:1588.8 3rd Qu.:0.8444 3rd Qu.:1043.3
## Max. :990.8 Max. :2409.1 Max. :0.9386 Max. :1984.4
##
## Harvest_index
## Min. :0.2569
## 1st Qu.:0.5143
## Median :0.5444
## Mean :0.5284
## 3rd Qu.:0.5847
## Max. :0.6348
## まずは、散布図を描いてデータをながめてみる。
colnames(df)## [1] "Plot" "Year"
## [3] "Rep" "Fertilizer"
## [5] "Variety" "GrowSta_N_stems"
## [7] "GrowSta_Plant_height" "GrowSta_SPAD"
## [9] "N_panicles" "N_grains_per_panicle"
## [11] "N_grains_total" "Ripening_rate"
## [13] "Thousand_grain_weight" "Grain_yield_gross"
## [15] "Grain_yield" "Sink_capacity"
## [17] "Sink_filling_rate" "Dry_weight_stems_leaves"
## [19] "Harvest_index"まずは、生育ステージのデータ（GrowSta_で始まる変数）だけで、確認してみる。
pairs(df[, 6:8],
col = df$Variety,
pch = as.numeric(df$Fertilizer),
oma = c(3, 3, 3, 15))
par(xpd = TRUE)
legend("topright",
legend = levels(df$Variety),
col = 1:nlevels(df$Variety),
pch = 1)
symbol_map <- c("circle", "square")
plot_ly(data = df,
x = ~GrowSta_N_stems,
y = ~GrowSta_Plant_height,
z = ~GrowSta_SPAD,
color = ~Variety,
symbol = ~Fertilizer,
symbols = symbol_map,
type = "scatter3d",
mode = "markers",
marker = list(size = 5))収穫時の形質のうち約半数を確認する。
pairs(df[, 9:14],
col = df$Variety,
pch = as.numeric(df$Fertilizer))
収穫時の形質のうち残りを確認する。
pairs(df[, 15:19],
col = df$Variety,
pch = as.numeric(df$Fertilizer))
収穫時の形質は次元が高すぎて把握が難しい。そこで主成分分析で要約してみる。計測尺度が異なるので必ず scale = TRUE とする。
res <- prcomp(df[, 6:19], scale = TRUE)主成分の寄与率を確認する。
plot(res, xlab = "PC")
寄与率で表す。
barplot(res$sdev^2 / sum(res$sdev^2),
xlab = "PC", ylab = "Contribution")
散布図を描く
pairs(res$x[,1:4],
col = df$Variety,
pch = as.numeric(df$Fertilizer))
biplotを描いて、主成分の意味を確認する。
biplot(res)
biplotでは点の形状や色を変えられず、その点の属性がわからない。そこで、ラベルの文字を変えてみる。
vf <- paste0(substr(df$Variety, 1, 1),
substr(df$Fertilizer, 1, 1))
biplot(res, xlabs = vf,
xlim = c(-0.5, 0.5),
ylim = c(-0.5, 0.5), cex = 0.8)
主成分ごとに箱ひげ図を描いてみる。
for(i in 1:4) {
boxplot(res$x[, i] ~ vf,
xlab = "Variety:Fertilizer",
ylab = paste0("PC", i))
}
品種効果・処理（施肥）効果・交互作用の確認（2023）
解析する形質を定める
colnames(df)## [1] "Plot" "Year"
## [3] "Rep" "Fertilizer"
## [5] "Variety" "GrowSta_N_stems"
## [7] "GrowSta_Plant_height" "GrowSta_SPAD"
## [9] "N_panicles" "N_grains_per_panicle"
## [11] "N_grains_total" "Ripening_rate"
## [13] "Thousand_grain_weight" "Grain_yield_gross"
## [15] "Grain_yield" "Sink_capacity"
## [17] "Sink_filling_rate" "Dry_weight_stems_leaves"
## [19] "Harvest_index"target <- "Grain_yield"
dt <- df %>% select(Year, Rep, Fertilizer, Variety) %>%
mutate(y = df[, target])
head(dt)## Year Rep Fertilizer Variety y
## 1 2023 1 Low_N Tachiayaka 713.7046
## 2 2023 1 Low_N Momiroman 544.6303
## 3 2023 1 Low_N Hokuriku193 990.7985
## 4 2023 1 Low_N Koshihikari 750.4786
## 5 2023 1 High_N Hokuriku193 795.3247
## 6 2023 1 High_N Momiroman 524.4486精籾収量について、品種ごとの平均を計算する。
VarMean <- tapply(dt$y, dt$Variety, mean)
VarMean## Hokuriku193 Koshihikari Momiroman Tachiayaka
## 821.8015 615.7702 664.2796 512.3007施肥水準ごとの平均を計算する。
FertMean <- tapply(dt$y, dt$Fertilizer, mean)
FertMean## High_N Low_N
## 661.1148 645.9612繰り返し（ブロック）ごとの平均を計算する。
RepMean <- tapply(dt$y, dt$Rep, mean)
RepMean## 1 2 3
## 622.4273 632.3664 705.8204品種と施肥水準の組合せごとの平均を計算する。
VarFertMean <- tapply(dt$y,
dt$Variety:dt$Fertilizer,
mean)
VarFertMean## Hokuriku193:High_N Hokuriku193:Low_N Koshihikari:High_N Koshihikari:Low_N
## 799.0363 844.5666 653.5172 578.0231
## Momiroman:High_N Momiroman:Low_N Tachiayaka:High_N Tachiayaka:Low_N
## 670.8348 657.7244 521.0709 503.5305品種と施肥水準を図示するために一度行列としてまとめる。
vf <- matrix(VarFertMean, nrow = 2)
vf## [,1] [,2] [,3] [,4]
## [1,] 799.0363 653.5172 670.8348 521.0709
## [2,] 844.5666 578.0231 657.7244 503.5305行と列に名前をつける。
rownames(vf) <- levels(dt$Fertilizer)
colnames(vf) <- levels(dt$Variety)
vf## Hokuriku193 Koshihikari Momiroman Tachiayaka
## High_N 799.0363 653.5172 670.8348 521.0709
## Low_N 844.5666 578.0231 657.7244 503.5305matplot関数を使って図示する。
matplot(vf, type = "b") # type = "b"は点と線の両方を描く
わかりやすいように軸のタイトル、ラベル、凡例をつける。
matplot(vf, type = "b",
pch = 1, ylim = c(0, max(vf)), xaxt = "n",
xlab = "Fertilizer",
ylab = target)
axis(1, at = 1:2, labels = rownames(vf))
legend("bottomright", legend = colnames(vf),
pch = 1, col = 1:4, cex = 1)
品種と施肥水準を逆にして棒グラフで図示してみる。
barplot(
vf,
beside = TRUE, # 並べて描画
col = c("skyblue", "orange"),# 棒の色
legend.text = TRUE, # 凡例を表示
args.legend = list(x = "bottomright"), # 凡例の位置
xlab = "Variety", # x軸ラベル
ylab = target # y軸ラベル
)
分散分析
分割区法分散分析のRでの実行は多少複雑である。一次因子の施肥水準と、二次因子の品種の効果、およびそれら交互作用を検定するが、一次因子の効果は一次誤差で、二次因子の効果と一次因子と二次因子の交互作用は二次誤差で検定される必要がある。
res <- aov(y ~
Fertilizer * Variety +
Error(Rep/Fertilizer),
data = dt)
summary(res)##
## Error: Rep
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 2 33196 16598
##
## Error: Rep:Fertilizer
## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizer 1 1378 1378 0.012 0.921
## Residuals 2 220587 110293
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Variety 3 298814 99605 4.459 0.0253 *
## Fertilizer:Variety 3 11000 3667 0.164 0.9185
## Residuals 12 268062 22339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ｒでは計算されないが、計算されている平均平方を用いて、反復間の差に対する検定も行うことができる。
# 反復のF値とｐ値
tmp <- summary(res)
(mss1 <- tmp$`Error: Rep`[[1]][1,3])## [1] 16598.24(mss2 <- tmp$`Error: Rep:Fertilizer`[[1]][2,3])## [1] 110293.4(df1 <- tmp$`Error: Rep`[[1]][1,1])## [1] 2(df2 <- tmp$`Error: Rep:Fertilizer`[[1]][2,1])## [1] 2fv_rep <- mss1/mss2; fv_rep## [1] 0.15049171 - pf(fv_rep, df1, df2)## [1] 0.86919362023年度の精籾収量の分散分析表は以下のようにまとめられる。
| 要因 | 自由度 | 平方和 | 平均平方 | F 値 | p 値 | |
|---|---|---|---|---|---|---|
| 反復 | 2 | 33196 | 16598 | 0.150 | 0.869 | |
| 施肥（1次） | 1 | 1378 | 1378 | 0.012 | 0.921 | |
| 1次誤差 | 2 | 220587 | 110293 | |||
| 品種（2次） | 3 | 298814 | 99605 | 4.459 | 0.0253 | * |
| 施肥・品種交互作用 | 3 | 11000 | 3667 | 0.164 | 0.919 | |
| 2次誤差 | 12 | 2680622 | 22339 | |||
2023年度の収量の分散分析では、一次因子である施肥効果は認められない。また、二次因子の品種効果は認められるが、施肥・品種の交互作用は認められない。これは、グラフをみると「北陸193号」以外は高窒素区で収量が上昇していたが、その効果は顕著でなかった。一方、品種間で大きな差があるのは明らかであった。交互作用は、4品種が高窒素区・低窒素区間で順位も変わらず、有意性は認められなかった。
主成分得点の分散分析
複数の形質から計算される主成分得点について、品種や処理の効果を含めた分散分析を行うこともできる。まず、データを準備する。
res <- prcomp(df[, 6:19], scale = TRUE)
target <- 1
dt <- df %>% select(Year, Rep, Fertilizer, Variety) %>%
mutate(y = res$x[,target])
head(dt)## Year Rep Fertilizer Variety y
## 1 2023 1 Low_N Tachiayaka 1.0183418
## 2 2023 1 Low_N Momiroman -2.8848191
## 3 2023 1 Low_N Hokuriku193 0.3037691
## 4 2023 1 Low_N Koshihikari -0.6277118
## 5 2023 1 High_N Hokuriku193 -1.1624640
## 6 2023 1 High_N Momiroman -2.0189045res <- aov(y ~
Fertilizer * Variety +
Error(Rep/Fertilizer),
data = dt)
summary(res)##
## Error: Rep
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 2 2.187 1.094
##
## Error: Rep:Fertilizer
## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizer 1 0.013 0.013 0.001 0.975
## Residuals 2 20.510 10.255
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Variety 3 83.73 27.911 13.340 0.000396 ***
## Fertilizer:Variety 3 0.69 0.232 0.111 0.952206
## Residuals 12 25.11 2.092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1年次効果を入れた分散分析
2023年度と2024年度では、施肥によって制御する対象が窒素からリンに変わった。また、2022年度と2023年度では、品種も異なる。そこで、ここでは2019年度から2022年度の圃場試験実習の結果をもとに年次効果を含めた分散分析について説明する。2019年度から2022年度では、年度ごとに配置は異なるが、同じ構造の試験であったので、年次間の比較が行うことができる。しかし、構造モデルはかなり複雑になる。
構造モデル
試験区のある形質の測定値を \(Y_{ijkl}\) とすると， \(Y_{ijkl}\) の統計（構造）モデルは，1 次因子水準数 \(a = 2\)，2 次因子水準数 \(b = 4\)，反復 \(r = 3\)，年次 \(c=4\) として，
\[ Y_{ijkl} = \mu + \alpha_i + \gamma_l + (\alpha \gamma)_{il} + ( \rho\gamma)_{kl} + e_{ikl}^1 + \beta_j + (\alpha \beta)_{ij} + (\beta\gamma)_{jl} + (\alpha \beta \gamma)_{ijl} + e_{ijkl}^2, \\ i=1,\ldots,a, j=1,\ldots,b, k=1,\ldots,r, l=1,\ldots,c \] とおける．ただし，\(\mu\) は総平均，\(\alpha_i\) は施肥効果，\(\gamma_l\) は年次効果，\((\alpha\gamma)_{il}\) は 施肥と年次の交互作用，\((\rho\gamma)_{kl}\) はブロック効果で、年次でブロックが異なるので交互作用の形で取り出した。\(\beta_j\) は品種効果，\((\alpha\beta)_{ij}\) は 施肥と品種の交互作用，\((\beta\gamma)_{jl}\) は 品種と年次の交互作用，\(e_{ikl}^1\) は 1 次誤差、\(e_{ijkl}^2\) は 2 次誤差、\((\alpha\beta\gamma)_{ijl}\) は 施肥と品種と年次の3次交互作用である。ブロック効果は年次をまたいで取り出すことはできないので、年次とブロックの交互作用として取り出し、\(\rho\gamma_{kl} \sim N(0, \sigma_B^2)\) を満たす変量（確率変数）とした。誤差は，\(e_{ikl}^1 \sim N(0, \sigma^2_1)\)、\(e^2_{ijkl} \sim N(0, \sigma^2_2)\) を満たす互いに独立な変量であるとした。
データの視覚化
2019年から2022年度のデータを読み込み、品種・施肥水準・年次の組合せごとに平均を求めて視覚化する。
re19.22 <- read.csv("riceexp19to22.csv")
re19.22 <- re19.22 %>% mutate(Plot = factor(Plot),
Year = factor(Year),
Rep = factor(Rep),
Fertilizer = factor(Fertilizer),
Variety = factor(Variety))
summary(re19.22)## Plot Year Rep Fertilizer Variety N_panicles
## 1 : 4 2019:24 1:32 High_N:48 Calotoc :24 Min. :124.4
## 2 : 4 2020:24 2:32 Low_N :48 Momiroman :24 1st Qu.:177.2
## 3 : 4 2021:24 3:32 Teqing :24 Median :200.0
## 4 : 4 2022:24 Yamadawara:24 Mean :207.8
## 5 : 4 3rd Qu.:233.3
## 6 : 4 Max. :304.4
## (Other):72
## Grain_yield Dry_weight_stems_leaves Harvest_index
## Min. : 78.35 Min. : 792.1 Min. :0.04213
## 1st Qu.:341.33 1st Qu.:1176.4 1st Qu.:0.24354
## Median :524.17 Median :1361.3 Median :0.35311
## Mean :504.61 Mean :1359.3 Mean :0.34169
## 3rd Qu.:664.38 3rd Qu.:1520.0 3rd Qu.:0.45123
## Max. :900.13 Max. :1848.7 Max. :0.62020
## target <- "Grain_yield"
dt <- re19.22 %>%
select(Year, Rep, Fertilizer, Variety) %>%
mutate(y = re19.22[, target])
head(dt)## Year Rep Fertilizer Variety y
## 1 2022 1 Low_N Calotoc 288.0794
## 2 2022 1 Low_N Yamadawara 595.0599
## 3 2022 1 Low_N Teqing 727.3062
## 4 2022 1 Low_N Momiroman 403.6212
## 5 2022 1 High_N Momiroman 522.1633
## 6 2022 1 High_N Calotoc 278.6429まずは、組合せごとの平均を計算する。
VarFertYearMean <- tapply(dt$y,
dt$Variety:dt$Fertilizer:dt$Year,
mean)
VarFertYearMean## Calotoc:High_N:2019 Calotoc:High_N:2020 Calotoc:High_N:2021
## 353.1543 379.2207 114.1943
## Calotoc:High_N:2022 Calotoc:Low_N:2019 Calotoc:Low_N:2020
## 256.4652 325.8597 286.2260
## Calotoc:Low_N:2021 Calotoc:Low_N:2022 Momiroman:High_N:2019
## 167.9413 279.3851 311.4990
## Momiroman:High_N:2020 Momiroman:High_N:2021 Momiroman:High_N:2022
## 508.2113 379.5093 571.3807
## Momiroman:Low_N:2019 Momiroman:Low_N:2020 Momiroman:Low_N:2021
## 353.3203 697.1863 575.9860
## Momiroman:Low_N:2022 Teqing:High_N:2019 Teqing:High_N:2020
## 448.8584 637.1487 737.8623
## Teqing:High_N:2021 Teqing:High_N:2022 Teqing:Low_N:2019
## 455.7623 825.4517 548.0263
## Teqing:Low_N:2020 Teqing:Low_N:2021 Teqing:Low_N:2022
## 767.2623 403.4193 761.7230
## Yamadawara:High_N:2019 Yamadawara:High_N:2020 Yamadawara:High_N:2021
## 658.2533 648.7797 454.7887
## Yamadawara:High_N:2022 Yamadawara:Low_N:2019 Yamadawara:Low_N:2020
## 704.8038 666.3393 667.6717
## Yamadawara:Low_N:2021 Yamadawara:Low_N:2022
## 566.1410 635.7732行列のかたちにする。
vfy <- matrix(VarFertYearMean, nrow = 4)
rownames(vfy) <- levels(dt$Year)
colnames(vfy) <- levels(dt$Variety:dt$Fertilizer)
vfy## Calotoc:High_N Calotoc:Low_N Momiroman:High_N Momiroman:Low_N
## 2019 353.1543 325.8597 311.4990 353.3203
## 2020 379.2207 286.2260 508.2113 697.1863
## 2021 114.1943 167.9413 379.5093 575.9860
## 2022 256.4652 279.3851 571.3807 448.8584
## Teqing:High_N Teqing:Low_N Yamadawara:High_N Yamadawara:Low_N
## 2019 637.1487 548.0263 658.2533 666.3393
## 2020 737.8623 767.2623 648.7797 667.6717
## 2021 455.7623 403.4193 454.7887 566.1410
## 2022 825.4517 761.7230 704.8038 635.7732matplot関数でグラフを描く。
matplot(vfy, type = "b")
軸タイトル、目盛り、凡例を加えて、再度描く。
matplot(vfy, type = "b", pch = 1:2,
col = rep(1:7, each = 2), lty = 1:2,
ylim = c(0, max(vfy)), xaxt = "n",
xlab = "Year", ylab = "Yield", )
axis(1, 1:nrow(vfy), labels = rownames(vfy))
legend("bottomright", legend = levels(dt$Variety),
lty = 1, col = 1:7, cex = 0.7)
legend("bottomleft", legend = levels(dt$Fertilizer),
lty = 1:2, pch = 1:2, cex = 0.7)
年次にもよるが、施肥水準による効果は、年次の効果や、品種の効果に比べると、あまり顕著ではない。
分散分析
ここからは、今年度のデータを除いて、分散分析を行ってみる。
年次効果も含めた分割区法分散分析を行う。
res <- aov(y ~ Year * Fertilizer * Variety
+ Error((Year + Year:Rep)/Fertilizer),
data = dt)
summary(res)##
## Error: Year
## Df Sum Sq Mean Sq
## Year 3 565467 188489
##
## Error: Year:Rep
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 8 91900 11487
##
## Error: Year:Rep:Fertilizer
## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizer 1 2242 2242 0.392 0.5486
## Year:Fertilizer 3 63329 21110 3.693 0.0619 .
## Residuals 8 45733 5717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Variety 3 2134493 711498 138.916 < 2e-16 ***
## Year:Variety 9 415933 46215 9.023 8.12e-08 ***
## Fertilizer:Variety 3 46689 15563 3.039 0.0379 *
## Year:Fertilizer:Variety 9 93357 10373 2.025 0.0567 .
## Residuals 48 245845 5122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ｒでは計算されないが、反復間の差や年次に対する検定も可能である。
tmp <- summary(res)
(mss1 <- tmp[[1]][[1]][1,3])## [1] 188489.1(mss2 <- tmp[[3]][[1]][3,3])## [1] 5716.584(df1 <- tmp[[1]][[1]][1,1])## [1] 3(df2 <- tmp[[3]][[1]][3,1])## [1] 8fv_year <- mss1/mss2; fv_year## [1] 32.972321 - pf(fv_year, df1, df2)## [1] 7.479304e-05(mss1 <- tmp[[2]][[1]][1,3])## [1] 11487.47(mss2 <- tmp[[3]][[1]][3,3])## [1] 5716.584(df1 <- tmp[[2]][[1]][1,1])## [1] 8(df2 <- tmp[[3]][[1]][3,1])## [1] 8fv_yearrep <- mss1/mss2; fv_yearrep## [1] 2.0094981 - pf(fv_yearrep, df1, df2)## [1] 0.17168432019～2022年度の収量の分散分析表は以下のようにまとめられる。
| 要因 | 自由度 | 平方和 | 平均平方 | F 値 | p 値 | |
|---|---|---|---|---|---|---|
| 年次 | 3 | 565467 | 188489 | 32.97 | 7e-15 | *** |
| 反復(年次) | 8 | 91900 | 11487 | 2.009 | 0.172 | |
| 施肥（1次） | 1 | 2242 | 2242 | 0.392 | 0.5486 | |
| 施肥・年次交互作用 | 3 | 63329 | 21110 | 3.693 | 0.0619 | . |
| 1次誤差 | 8 | 45733 | 5717 | |||
| 品種（2次） | 3 | 2134493 | 711498 | 138.916 | < 2e-16 | *** |
| 品種・年次交互作用 | 9 | 415933 | 46215 | 9.023 | 8.12e-08 | *** |
| 施肥・品種交互作用 | 3 | 46689 | 15563 | 3.039 | 0.0379 | * |
| 施肥・品種・年次3次交互作用 | 9 | 93357 | 10373 | 2.025 | 0.0567 | . |
| 2次誤差 | 48 | 245845 | 5122 | |||
施肥効果は認められなかったが、品種効果、品種・年次交互作用、施肥・品種交互作用が認められた。3次交互作用は微妙だが認められなかった。交互作用が認められた組のグラフを描く。
まず、年次と品種の組合せごとに平均を求める。
VarYear <- tapply(dt$y,
dt$Variety:dt$Year,
mean)
head(VarYear)## Calotoc:2019 Calotoc:2020 Calotoc:2021 Calotoc:2022 Momiroman:2019
## 339.5070 332.7233 141.0678 267.9251 332.4097
## Momiroman:2020
## 602.6988行列のかたちにする。
vy <- matrix(VarYear, nrow = 4)
rownames(vy) <- levels(dt$Year)
colnames(vy) <- levels(dt$Variety)
vy## Calotoc Momiroman Teqing Yamadawara
## 2019 339.5070 332.4097 592.5875 662.2963
## 2020 332.7233 602.6988 752.5623 658.2257
## 2021 141.0678 477.7477 429.5908 510.4648
## 2022 267.9251 510.1195 793.5874 670.2885図示する。
matplot(vy, type = "b", pch = 1, lty = 1, col = 1:4,
ylim = c(0, max(vy)), xaxt = "n",
xlab = "Year", ylab = "Yield", )
axis(1, 1:4, labels = rownames(vy))
legend("bottomleft", legend = colnames(vy),
pch = 1, col = 1:4, cex = 0.8)
品種と施肥水準の組合せごとに平均を求める。
VarFertMean <- tapply(dt$y,
dt$Variety:dt$Fertilizer,
mean)
head(VarFertMean)## Calotoc:High_N Calotoc:Low_N Momiroman:High_N Momiroman:Low_N
## 275.7586 264.8530 442.6501 518.8378
## Teqing:High_N Teqing:Low_N
## 664.0563 620.1077行列の形にする。
vf <- matrix(VarFertMean, nrow = 2)
rownames(vf) <- levels(dt$Fertilizer)
colnames(vf) <- levels(dt$Variety)
vf## Calotoc Momiroman Teqing Yamadawara
## High_N 275.7586 442.6501 664.0563 616.6564
## Low_N 264.8530 518.8378 620.1077 633.9813図示する。
matplot(vf, type = "b", pch = 1, lty = 1, col = 1:4,
ylim = c(0, max(vf)), xaxt = "n",
xlab = "Fertilizer", ylab = "Yield", )
axis(1, 1:2, labels = rownames(vf))
legend("bottomleft", legend = colnames(vf),
pch = 1, col = 1:4, cex = 0.7)
重回帰分析による関連の解析（2023）
重回帰を用いて、生育期の形質と成熟期の形質間にどのような関係があるかを解析してみる。
重回帰分析を用いて、精籾収量を生育期の形質に回帰するモデルを作成する。
model <- lm(Grain_yield ~
GrowSta_N_stems + GrowSta_Plant_height + GrowSta_SPAD,
data = df)
summary(model)##
## Call:
## lm(formula = Grain_yield ~ GrowSta_N_stems + GrowSta_Plant_height +
## GrowSta_SPAD, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -247.42 -143.33 23.64 84.56 251.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 719.0750 375.8505 1.913 0.07015 .
## GrowSta_N_stems -0.1791 0.3774 -0.475 0.64015
## GrowSta_Plant_height -17.7069 5.8629 -3.020 0.00676 **
## GrowSta_SPAD 17.8323 10.5796 1.686 0.10742
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 155.3 on 20 degrees of freedom
## Multiple R-squared: 0.4208, Adjusted R-squared: 0.3339
## F-statistic: 4.843 on 3 and 20 DF, p-value: 0.01081草丈の効果が高度に有意であることがわかる。
しかし、用いられた品種には草丈の違いがあり、4品種がたまたまもつ草丈の違いが収量と関連しただけかもしれない。そこで、先に行った分散分析のモデルと重回帰モデルを融合して、草丈の有意性を確認することにする。
まず、先に行った分散分析と同じ構造をもつ回帰モデルを構築する。一次因子の誤差は、ここでは、施肥水準と反復の交互作用（Fertilizer:Rep）である。実際には、施肥水準と反復の主効果の検定は、施肥水準と反復の交互作用を誤差として行わなければならないことに注意する。
model <- lm(Grain_yield ~
Fertilizer + Variety + Rep +
Rep:Fertilizer + Fertilizer:Variety,
data = df)
anova(model)## Analysis of Variance Table
##
## Response: Grain_yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizer 1 1378 1378 0.0617 0.80806
## Variety 3 298814 99605 4.4589 0.02526 *
## Rep 2 33196 16598 0.7430 0.49634
## Fertilizer:Rep 2 220587 110293 4.9374 0.02725 *
## Fertilizer:Variety 3 11000 3667 0.1641 0.91846
## Residuals 12 268062 22339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1上のモデルに、先の重回帰モデルを融合する。
model <- lm(Grain_yield ~
Fertilizer + Variety + Rep +
Rep:Fertilizer + Fertilizer:Variety +
GrowSta_N_stems + GrowSta_Plant_height + GrowSta_SPAD,
data = df)
anova(model)## Analysis of Variance Table
##
## Response: Grain_yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizer 1 1378 1378 0.0718 0.79472
## Variety 3 298814 99605 5.1935 0.02353 *
## Rep 2 33196 16598 0.8654 0.45314
## GrowSta_N_stems 1 138113 138113 7.2013 0.02506 *
## GrowSta_Plant_height 1 48816 48816 2.5453 0.14509
## GrowSta_SPAD 1 2711 2711 0.1413 0.71566
## Fertilizer:Rep 2 118607 59303 3.0921 0.09502 .
## Fertilizer:Variety 3 18793 6264 0.3266 0.80631
## Residuals 9 172609 19179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1生育期の草丈は有意でなく、先ほどのモデルでは有意でなかった生育期の茎数が5%水準で有意であることがわかる。
効果の向きを確認する。
coef(model)## (Intercept) FertilizerLow_N
## 829.1611277 211.8852954
## VarietyKoshihikari VarietyMomiroman
## -29.3960706 -161.1714698
## VarietyTachiayaka Rep2
## -100.0674499 217.9012924
## Rep3 GrowSta_N_stems
## 134.5907079 -0.3426436
## GrowSta_Plant_height GrowSta_SPAD
## -16.9998249 15.1060247
## FertilizerLow_N:Rep2 FertilizerLow_N:Rep3
## -363.4028493 -139.2355250
## FertilizerLow_N:VarietyKoshihikari FertilizerLow_N:VarietyMomiroman
## -168.3527995 -93.2544689
## FertilizerLow_N:VarietyTachiayaka
## -72.2696723生育期の茎数の回帰係数は負であり、茎数が多くなるほど、精籾収量は低下すると考えられる。