# ーーーーーー【偏差分割のグラフ化】ーーーー−ーーー
#
#  データのもつ全偏差を処理偏差と誤差偏差に分割し,
#  それぞれの偏差をインデックス・プロットの図上に
#  図示するRスクリプト.
# ーーーーーーーーーーーーーーーーーーーーーーーーー

#−−−−【偏差分割のグラフ化】−−−−

results <- read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/yields.txt", header =T)

# Michael J. Crawley (2007), The R Book のテストデータ
# http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/yields.txt

attach(results)
data <- c(sand, clay, loam)

# インデックス・プロット
index <- 1:30
plot(data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 全偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
for (i in 1:30) lines(c(i,i), c(data[i], mean(data)), lty=2, lwd=2)
abline(h=mean(data), lwd=4)
text(25, 5, "total deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 処理効果を仮定しない「帰無仮説」は,データのばらつきを
# この「全偏差の図示」にしたがって説明しようと宣言する.

# 処理偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
abline(h=mean(data), lwd=4)
lines(c(1,10), c(mean(sand),mean(sand)), lwd=2)
lines(c(11,20), c(mean(clay),mean(clay)), lwd=2)
lines(c(21,30), c(mean(loam),mean(loam)), lwd=2)
for (i in 1:10) lines(c(i,i), c(mean(sand), mean(data)), lty=2)
for (i in 11:20) lines(c(i,i), c(mean(clay), mean(data)), lty=2)
for (i in 21:30) lines(c(i,i), c(mean(loam), mean(data)), lty=2)
text(23, 5, "treatment deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 誤差偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
lines(c(1,10), c(mean(sand),mean(sand)), lwd=2)
lines(c(11,20), c(mean(clay),mean(clay)), lwd=2)
lines(c(21,30), c(mean(loam),mean(loam)), lwd=2)
for (i in 1:10) lines(c(i,i), c(data[i], mean(sand)), lty=3)
for (i in 11:20) lines(c(i,i), c(data[i], mean(clay)), lty=3)
for (i in 21:30) lines(c(i,i), c(data[i], mean(loam)), lty=3)
text(25, 5, "error deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 処理効果をモデルに含む「対立仮説」は,データのばらつきを
# この「誤差偏差の図示」にしたがって説明しようと宣言する.
# まとめて描画(2×2)

par(mfrow=c(2,2))

# インデックス・プロット
index <- 1:30
plot(data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 全偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
for (i in 1:30) lines(c(i,i), c(data[i], mean(data)), lty=2, lwd=2)
abline(h=mean(data), lwd=4)
text(25, 5, "total deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 処理偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
abline(h=mean(data), lwd=4)
lines(c(1,10), c(mean(sand),mean(sand)), lwd=2)
lines(c(11,20), c(mean(clay),mean(clay)), lwd=2)
lines(c(21,30), c(mean(loam),mean(loam)), lwd=2)
for (i in 1:10) lines(c(i,i), c(mean(sand), mean(data)), lty=2)
for (i in 11:20) lines(c(i,i), c(mean(clay), mean(data)), lty=2)
for (i in 21:30) lines(c(i,i), c(mean(loam), mean(data)), lty=2)
text(23, 5, "treatment deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")

# 誤差偏差の図示
plot(index, data, ylim=c(3,18), pch=(15+(index>10)+(index>20)))
lines(c(1,10), c(mean(sand),mean(sand)), lwd=2)
lines(c(11,20), c(mean(clay),mean(clay)), lwd=2)
lines(c(21,30), c(mean(loam),mean(loam)), lwd=2)
for (i in 1:10) lines(c(i,i), c(data[i], mean(sand)), lty=3)
for (i in 11:20) lines(c(i,i), c(data[i], mean(clay)), lty=3)
for (i in 21:30) lines(c(i,i), c(data[i], mean(loam)), lty=3)
text(25, 5, "error deviation")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")