#################################################
#豊田秀樹(1994) 違いを見ぬく統計学　実験計画と分散分析入門　講談社
#第8章 p.222-246
#################################################

#---------------------------#
#データの読み込み
#---------------------------#
hotel5 <- read.csv("http://blue.zero.jp/yokumura/R/ancova/hotel5.csv",header=T)
hotel5
ancova <- data.frame(
	subset(hotel5,季節=="夏",select=c("料金","都市")),
	subset(hotel5,季節=="冬",select=c("料金","都市"))
	)
names(ancova) <- c("夏料金","都市","冬料金","都市")
ancova <- ancova[,1:3]
ancova

#---------------------------#
#図 8-5 (p.234)を描画
#---------------------------#
plot(夏料金~冬料金, data=ancova, type="n")
lines(夏料金~冬料金, data=ancova,type="p", col="red", subset=都市=="シカゴ")
lines(夏料金~冬料金, data=ancova,type="p", col="blue", subset=都市=="オーランド")
abline(lm(夏料金~冬料金, data=ancova))

graphics.off()

library(car)		#scatterplot() を利用すると簡単に散布図が描ける
scatterplot(夏料金~冬料金|都市,data=ancova)
graphics.off()


##夏料金は，シカゴとオーランドで差があるか？
summary(aov(夏料金~都市,data=ancova))

##Step1: 補助因子と従属変数の相関を確認
cor.test(~夏料金+冬料金,data=ancova)
cor.test(~夏料金+冬料金,data=ancova,subset=都市=="シカゴ")
cor.test(~夏料金+冬料金,data=ancova,subset=都市=="オーランド")


##Step2: 回帰係数は群間で等しいか？
mod1 <- lm(夏料金 ~ 都市 * 冬料金,data=ancova)
mod2 <- lm(夏料金 ~ 都市 + 冬料金,data=ancova)
mod3 <- lm(夏料金 ~ 冬料金,data=ancova)
anova(mod2,mod1)		#図表 8-6のtypeIIIの平方和


##Step3: 切片は群間で等しいか？
library(gmodels)		#ci()のため R version 2.0 未満の場合，library(gregmisc)
anova(mod3,mod2)		#図表 8-7のtypeIIIの平方和
ci(mod2)			#冬料金の回帰係数の推定値は1.27，標準誤差は.29


##Step4: 切片の推定値
ci(mod2)	#オーランド
		#都市 が オーランドの場合
		#夏料金 = -44.407813（切片） + 0 (都市)  +  1.271769 * 冬料金
		#
		#都市 が シカゴの場合
		#夏料金 = -44.407813（切片） +  44.451891 (都市)  +  1.271769 * 冬料金
		#夏料金 =  0.044078  (切片)  +  1.271769 * 冬料金
		


#reference を変更して確認する
ancova$都市	#levelsは，オーランド，シカゴの順になっている
		#levelsで第1番目に標示されるものが，referenceになるので，
		#referenceを変更することでも
		#各都市の切片の推定値を算出することが出来る．
ancova$都市 <- relevel(ancova$都市,ref="シカゴ")
ancova$都市

mod4 <- lm(夏料金 ~ 都市 + 冬料金,data=ancova)
ci(mod4)	#シカゴ



##共分散分析の結果を図示する
Chicago <- subset(ancova,都市=="シカゴ")
Orland  <- subset(ancova,都市=="オーランド")

plot(夏料金~冬料金,data=ancova,type="n")	#xlimなどを自動で指定する
lines(夏料金~冬料金,data=Chicago,type="p",col="red")
lines(夏料金~冬料金,data=Orland,type="p",col="blue")

abline(-44.407813, 1.271769, col="blue")
abline(-44.407813 + 44.451891, 1.271769, col="red")

graphics.off()