# ーーーーーー【偏差分割のグラフ化】ーーーー−ーーー
#
# データのもつ全偏差を処理偏差と誤差偏差に分割し，
# それぞれの偏差をインデックス・プロットの図上に
# 図示するＲスクリプト．
# ーーーーーーーーーーーーーーーーーーーーーーーーー
#−−−−【偏差分割のグラフ化】−−−−
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")
par(mfrow=c(1,1))
# まとめて描画（1×2）
par(mfrow=c(1,2))
# 帰無仮説の世界
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, "H0: null")
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, "H1: alternative")
text(5, 18, "sand")
text(15, 18, "clay")
text(25, 18, "loam")
par(mfrow=c(1,1))