#################################################
#古川俊之 監修　丹後俊郎著　(1993)
#新版 医学への統計学 朝倉書店
#第6章
#第13章
#################################################

#-----------------------------------------------#
#6.2.3 2つの母回帰直線の差の検定 p.105-112
#-----------------------------------------------#

#----------------------------#
#Step1: データの読み込み
#----------------------------#
library(gmodels)		#ci()のため R version 2.0 未満の場合，library(gregmisc)
library(car)			#scatterplot() を利用すると簡単に散布図が描ける

dat1 <- read.csv("http://blue.zero.jp/yokumura/R/ancova/tango_1.csv",header=T)
dat1
dat1$地区 <- relevel(dat1$地区, ref="千葉県")	#千葉を第1群とする
attach(dat1)


#----------------------------#
#Step2: 準備
#----------------------------#

##(a) Xの平均
(Xbar <- tapply(NO2,地区,mean))
(Xbart <- mean(NO2))


##(b) Yの平均
(Ybar <- tapply(有症者率,地区,mean))
(Ybart<-mean(有症者率))


##(c) XYの積和
(SSxy <- tapply(NO2 * 有症者率, 地区, sum) - 
	tapply(NO2, 地区, sum) * tapply(有症者率, 地区, sum)/tapply(地区, 地区, length))
(SSxyt <- sum(NO2 * 有症者率) - sum(NO2) * sum(有症者率) / nrow(dat1))


##(d) Xの平方和
(SSx <- tapply(NO2^2, 地区, sum) - tapply(NO2, 地区, sum)^2 /  tapply(地区, 地区, length))
(SSxt <- sum(NO2^2) - sum(NO2)^2 / nrow(dat1))

##(e) Yの平方和
(SSy <- tapply(有症者率^2, 地区, sum) - tapply(有症者率, 地区, sum)^2 /  tapply(有症者率, 地区, length))
(SSyt <- sum(有症者率^2) - sum(有症者率)^2 / nrow(dat1))


#----------------------------#
#Step3: 回帰分析（2群を別々に）
#----------------------------#

##(a) 回帰係数（勾配係数）
(slope <- SSxy/SSx)


##(b) 切片
Ybar - slope * Xbar


##(c) 相関係数
SSxy / sqrt(SSx * SSy)



#----------------------------#
#Step4: 回帰直線が平行か否かの検定（6.27, 28, 29式）
#----------------------------#
(delta1 <- sum( SSy - SSxy^2/SSx ))
(delta2 <- sum(SSy) - sum(SSxy)^2 / sum(SSx))
(Fvalue <- (delta2-delta1) / (delta1/(nrow(dat1)-4)))
qf(p=.05, df1=1, df2=nrow(dat1)-4, lower.tail=F)



#----------------------------#
#Step5: 切片の相当性の検定（6.36, 37, 38式）
#----------------------------#
(delta3 <- SSyt-(SSxyt)^2/SSxt)
(delta3-delta2) / (delta2/(nrow(dat1)-3))
qf(p=.01, df1=1, df2=nrow(dat1)-3, lower.tail=F)



#----------------------------#
#Step6: 回帰係数の推定 (6.35, 39, 40式)
#----------------------------#
(bE <- sum(SSxy)/sum(SSx))
(intercept <- Ybar - bE * Xbar)



#----------------------------#
#Step7: 信頼限界 (6.43式)
#----------------------------#
(lower <- (intercept[[2]] - intercept[[1]]) - 
	qt(p=(.05/2), df=(nrow(dat1)-3), lower.tail=F) * 
	sqrt(
		(sum(1/tapply(地区,地区,length)) + (Xbar[[1]]-Xbar[[2]])^2/(sum(SSx))) *  
		delta2 / (nrow(dat1)-3)
	))

(upper <- (intercept[[2]] - intercept[[1]]) +
	qt(p=(.05/2), df=(nrow(dat1)-3), lower.tail=F) * 
	sqrt(
		(sum(1/tapply(地区,地区,length)) + (Xbar[[1]]-Xbar[[2]])^2/(sum(SSx))) *  
		delta2 / (nrow(dat1)-3)
	))




#----------------------------#
#Step8: 関数使用
#----------------------------#
scatterplot(有症者率 ~ NO2 | 地区,data=dat1)

cor.test(~ NO2 + 有症者率)
cor.test(~ NO2 + 有症者率,subset=地区=="千葉県")
cor.test(~ NO2 + 有症者率,subset=地区=="岡山県")


reg1 <- lm(有症者率 ~ NO2,data=dat1)
reg2 <- lm(有症者率 ~ 地区 + NO2,data=dat1)
reg3 <- lm(有症者率 ~ 地区 * NO2,data=dat1)


##回帰直線が平行か否かの検定
anova(reg2,reg3)


##切片の相当性の検定
anova(reg1,reg2)


##回帰係数の推定
anova(reg2)
summary(reg2)
ci(reg2)
slope.chiba <- ci(reg2)[1,1] + 0
slope.okayama <- ci(reg2)[1,1] + ci(reg2)[2,1]

plot(有症者率 ~ NO2,pch=as.numeric(地区),data=dat1)
abline(slope.chiba,ci(reg2)[3,1]  ,col="red")
abline(slope.okayama,ci(reg2)[3,1],col="blue")


graphics.off()
detach(dat1)




#-----------------------------------------------#
#13.1 共分散分析 p.240-244
#-----------------------------------------------#


##(a) データの読み込み
dat2 <- read.csv("http://blue.zero.jp/yokumura/R/ancova/tango_2.csv",header=T)


##(b) 地区ごとの血圧の差
boxplot(血圧~地区,data=dat2)


##(b) 血圧の比較 表 107下(p.241)
attach(dat2)
tapply(地区,地区,length)
tapply(血圧,地区,mean)
tapply(血圧,地区,sd)
detach(dat2)

t.test(血圧~地区,data=dat2, var.equal=T)


##(c) 2地区における血圧と年齢の関連 図82(p.242)
scatterplot(血圧~年齢|地区, data=dat2)


##(d) 年齢の比較 表
attach(dat2)
tapply(地区,地区,length)
tapply(年齢,地区,mean)
tapply(年齢,地区,sd)
detach(dat2)

t.test(年齢~地区,data=dat2, var.equal=T)


##(e) 共分散分析
cor.test(~血圧+年齢, data=dat2)
cor.test(~血圧+年齢, subset=地区=="A", data=dat2)
cor.test(~血圧+年齢, subset=地区=="B", data=dat2)

reg1 <- lm(血圧~年齢, data=dat2)
reg2 <- lm(血圧~年齢 + 地区, data=dat2)
reg3 <- lm(血圧~年齢 * 地区, data=dat2)

anova(reg2,reg3)		#回帰係数が平行か否か？
anova(reg1,reg2)		#切片が等しいか否か？

anova(reg1)
ci(reg1)


graphics.off()