#----------------------------------------------------
#永田（2003）第8章 対応がある場合の母平均の差の検定
#----------------------------------------------------


######################################
#対応がある場合の母平均の検定（検出力）
#(1)u=/ u_0　（精確値）
######################################

power.t1 <- function(sig.level,delta,n){
	pt(qt(p=sig.level/2,df=n-1,lower.tail=T),df=n-1,ncp=delta*sqrt(n),lower.tail=T)+
	pt(qt(p=sig.level/2,df=n-1,lower.tail=F),df=n-1,ncp=delta*sqrt(n),lower.tail=F)
	}



######################################
#対応がある場合の母平均の検定（検出力）
#(2)u > u_0　（精確値）
######################################

power.t2 <- function(sig.level,delta,n){
	pt(qt(p=sig.level,df=n-1,lower.tail=F),df=n-1,ncp=delta*sqrt(n),lower.tail=F)
	}



######################################
#対応がある場合の検定（検出力）
#(3)u < u_0　（精確値）
######################################
power.t3 <- function(sig.level,delta,n){
	pt(qt(p=sig.level,df=n-1,lower.tail=T),df=n-1,ncp=delta*sqrt(n),lower.tail=T)
	}



######################################
#対応がある場合の検定（サンプルサイズ）
#(1)u =/ u_0
######################################
samplesize.t1 <- function(power,delta, sig.level){
	((qnorm(sig.level/2,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level/2,lower.tail=F)^2)/2
}

samplesize.t1(power=.9,sig.level=.05,delta=1)


######################################
#対応がある場合の母平均の検定（サンプルサイズ）
#(2)u > u_0
######################################
samplesize.t2 <- function(power,delta, sig.level){
	((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level,lower.tail=F)^2)/2
}
samplesize.t2(power=.95,sig.level=.05,delta=.5)


######################################
#対応がある場合の母平均の検定（サンプルサイズ）
#(3)u < u_0　（精確値）（2）と同一
######################################
samplesize.t3 <- function(power,delta, sig.level){
	((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level,lower.tail=F)^2)/2
}
samplesize.t3(power=.8,sig.level=.05,delta=-.5)




#ex 8.1
samplesize.t1(power=.9,sig.level=.05,delta=1)
power.t1(n=13,sig.level=.05,delta=1)
power.t1(n=12,sig.level=.05,delta=1)


#ex 8.2
samplesize.t2(power=.95,sig.level=.05,delta=0.5)
power.t2(n=45,sig.level=.05,delta=0.5)
power.t2(n=44,sig.level=.05,delta=0.5)


##ex 8.3
#ex 8.2
samplesize.t3(power=.80,sig.level=.05,delta=-0.5)
power.t3(n=27,sig.level=.05,delta=-0.5)
power.t3(n=26,sig.level=.05,delta=-0.5)



##pr 8.1
x1 <- c(6.2,4.8,7.3,5.5,6.5,4.9,6.8,7.9,6.6,7.3)
x2 <- c(5.3,6.2,5.9,7.3,8.4,7.3,6.9,7.6,8.5,8.1)
t.test(x1,x2,paired=T)


##pr 8.2
x1 <- c(10.8,11.2,9.7,9.9,12.0,9.6,10.5,10.7,10.1)
x2 <- c(10.2,10.1,9.9,8.2,10.2,9.4,10.4,10.0,10.5)
t.test(x1,x2,paired=T,alternative="greater")

##pr 8.3
x1 <- c(19,19,21,19,22,19,20,21)
x2 <- c(21,22,16,22,25,18,24,23)
t.test(x1,x2,paired=T,alternative="less")

##pr 8.4
power.t1(sig.level=0.05,delta=0.5,n=10)
power.t1(sig.level=0.05,delta=0.5,n=20)

##pr 8.5
power.t2(sig.level=0.05,delta=1,n=9)
power.t2(sig.level=0.05,delta=1.5,n=9)

##pr 8.6
power.t3(sig.level=0.05,delta=-1,n=8)
power.t3(sig.level=0.05,delta=-1,n=16)