4年間にわたって行った データ表 を用いると， 気象条件などの年ごとの違いによる効果がわかる． 年次効果を γ_k で表すとすると，

「日本晴」収量の年次効果と施肥水準によるプロット（1996 〜 2000年度）

　上の分散分析表をみると，年次効果と施肥量効果は主効果と交互作用ともに有意 となったが，密度に関しては有意とはならなかった．どのような効果が存在するかをみるため，年次と施肥量水準の 組み合わせごとの平均を図にプロットした．
　これをみると， 年次ごとに収量に変動があることがわかる．また， 無肥料区(fert1)と少肥料区(fert2)は年次ごとの変動が似ていて，少肥料区 の方が収量が高くなっているが，多肥料区(fert3)は2000年のように高収量 のときもあるが，1996年や1997年のようにそれほど高い収量が得られない場合 もあり，これが交互作用として検出されている．これは，多肥料区で繁茂しすぎて 倒伏が起こった年があったためと考えられる．

R のコマンド

• 4 年次収量データ分散分析

  yield <- rice.dfr$gy	# gy 列を取り出す
  dens <- factor(rice.dfr$density)	#密度水準ラベル化
  fert <- factor(rice.dfr$fert)
  year <- factor(rice.dfr$year)
  rya.aov <- aov(yield ~ dens + fert + year +
  dens:fert + fert:year +	dens:year + dens:fert:year)
  	# 3 次交互作用まで計算
  summary(rya.aov)

• 「日本晴」収量の年次効果と施肥水準によるプロット
  x <- tapply(rice.dfr$gy, fert:year, mean)
  plot(1:4, x[1:4], type="b", lwd=2, cex=1.5, xaxt="n", xlab="year", ylab="yield",
         ylim=c(300,600), pch=0, lty=1,cex.lab=1.2, cex.axis=1.2, col="red")
  axis(1, 1:4, labels=levels(year), cex.axis=1.2)
  box(lwd=3)
  par(new=T)
  plot(1:4, x[5:8], type="b",lwd=2, cex=1.5,
         axes=F, xlab="", ylab="" , ylim=c(300,600), pch=2,lty=1, col="blue")
  par(new=T)
  plot(1:4, x[9:12], type="b", lwd=2, cex=1.5, axes=F, xlab="", ylab="" ,
         ylim=c(300,600), pch=3,lty=1, col="green")
  legend(locator(1), legend=c("fert1","fert2","fert3"), pch=c(0,2,3),
         col=c("red","blue","green"), cex=1.2)