#----------------------------------------------------
#永田（2003）第7章 2つの母平均の差の検定
#----------------------------------------------------

###############################
#2つの母平均の差の検定（検出力）
#(1)H_1 /= H_2（精確値）
###############################
power.t <- function(sig.level,delta,n1,n2){
	pt(qt(p=sig.level/2,df=n1+n2-2,lower.tail=T),
	df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=T)+
	pt(qt(p=sig.level/2,df=n1+n2-2,lower.tail=F),df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=F)
	}
power.t(sig.level=.05,delta=.6,n1=10,n2=8)


###############################
#2つの母平均の差の検定（検出力）
#(2)H_1 > H_2（精確値）
###############################
power.t2 <- function(sig.level,delta,n1,n2){
	pt(qt(p=sig.level,df=n1+n2-2,lower.tail=F),df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=F)
	}
power.t2(sig.level=.05,delta=.6,n1=10,n2=8)


###############################
#2つの母平均の差の検定（検出力）
#(3)H_1 < H_2（精確値）
###############################
power.t3 <- function(sig.level,delta,n1,n2){
	pt(qt(p=sig.level,df=n1+n2-2,lower.tail=T),
	df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=T)
	}
power.t3(sig.level=.05,delta=-.6,n1=10,n2=8)


###############################
#2つの母平均の差の検定（サンプルサイズ）
#(1)H_1 =/ H_2（精確値）
###############################
samplesize.t1 <- function(power,delta, sig.level){
	2*((qnorm(sig.level/2,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level/2,lower.tail=F)^2)/4
}

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


###############################
#2つの母平均の差の検定（サンプルサイズ）
#(2)H_1 > H_2（精確値）
###############################
samplesize.t2 <- function(power,delta, sig.level){
	2*((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level,lower.tail=F)^2)/4
}

samplesize.t2(power=.95,sig.level=.05,delta=1.5)


###############################
#2つの母平均の差の検定（サンプルサイズ）
#(3)H_1 < H_2（精確値）
###############################
samplesize.t3 <- function(power,delta, sig.level){
	2*((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	delta)^2+
	(qnorm(sig.level,lower.tail=F)^2)/4
}

samplesize.t3(power=.8,sig.level=.05,delta=-1.2)



##ex 7.2
delta <- seq(-1,1,by=.01)
plot(power.t(sig.level=.05,delta=delta,n1=10,n2=8)~delta,
type="l",ylim=c(0,1))
lines(power.t(sig.level=.05,delta=delta,n1=21,n2=16)~delta)
lines(power.t(sig.level=.05,delta=delta,n1=31,n2=25)~delta)

delta <- c(seq(-1,-.2,by=.2),seq(-.1,.1,by=.05),seq(.2,1,by=.2))
power.t(sig.level=.05,delta=delta,n1=10,n2=8)


##ex 7.4
delta <- seq(-1,1,by=.01)
plot(power.t2(sig.level=.05,delta=delta,n1=10,n2=8)~delta,
type="l",ylim=c(0,1))
lines(power.t2(sig.level=.05,delta=delta,n1=21,n2=16)~delta)
lines(power.t2(sig.level=.05,delta=delta,n1=31,n2=25)~delta)

delta <- c(seq(-1,-.2,by=.2),seq(-.1,.1,by=.05),seq(.2,1,by=.2))
power.t2(sig.level=.05,delta=delta,n1=10,n2=8)

##ex 7.6
delta <- seq(-1,1,by=.01)
plot(power.t3(sig.level=.05,delta=delta,n1=10,n2=8)~delta,
type="l",ylim=c(0,1))
lines(power.t3(sig.level=.05,delta=delta,n1=21,n2=16)~delta)
lines(power.t3(sig.level=.05,delta=delta,n1=31,n2=25)~delta)

delta <- c(seq(-1,-.2,by=.2),seq(-.1,.1,by=.05),seq(.2,1,by=.2))
power.t3(sig.level=.05,delta=delta,n1=10,n2=8)


##ex 7.7
samplesize.t1(power=.9,sig.level=.05,delta=1)
power.t(sig.level=.05,delta=1,n1=22,n2=22)


##ex 7.8
samplesize.t2(power=.95,sig.level=.05,delta=1.5)
power.t2(sig.level=.05,delta=1.5,n1=11,n2=11)
power.t2(sig.level=.05,delta=1.5,n1=10,n2=10)


##ex 7.9
samplesize.t3(power=.8,sig.level=.05,delta=-1.2)
power.t3(n1=10,n2=10,sig.level=.05,delta=-1.2)
power.t3(n1=9,n2=9,sig.level=.05,delta=-1.2)