################################################
#南風原（2002）第5章
################################################
#Step1：最尤法 p.129-31
################################################

##(a)尤度関数を定義 （5.3）式
ML <- function(PI){
	PI^3*((1-PI)^2)
	}


##(b)図5-1を描画
x <- seq(0,1,by=0.01)	#x軸を作成
mlx <- ML(x)		#y軸を作成
plot(mlx~x, type="l",
	xlab="π", ylab="尤度")



################################################
#Step2：相関係数の標準誤差 p.132-33
################################################

##(a)相関係数(r)と標本サイズ(N)より「標準誤差」を求める関数の定義 （5.7式）
r.sem <- function(r,N){
		(1-r^2)/N
		}

r.sem(.5,100)		#r=0.05 N=100の場合



################################################
#Step3：標準誤差に基づくサンプルサイズの決定 p.133-34
################################################

##(a)相関係数(r)と標準誤差(sem)より「サンプルサイズ」を求める関数の定義
power.cor <- function(r,sem){
		((1-r^2)^2)*(1/sem)^2
		}

power.cor(.5,.05)	#r=0.5 sem=0.05の場合
power.cor(0,.05)	#r=0   sem=0.05の場合



################################################
#Step4：検定の考え方 p.134-142
################################################

##(a)付表3 相関係数の棄却の限界値を再現 p.378
library(SuppDists)			##SuppDistsというライブラリを呼び出す
qPearson(p=.025, N=20, lower.tail=F)	#qPearson(p=有意水準(片側),N=サンプルサイズ(4以上),lower.tail=F 上側確率)
N <- c(
	seq(4,20,by=1), seq(22,40,by=2), seq(50,100,by=10), seq(200,400,by=100)
	)				#サンプルサイズを指定
	
x1 <- qPearson(p=.05,  N=N, lower.tail=F)	#有意水準.10（両側）
x2 <- qPearson(p=.025, N=N, lower.tail=F)	#有意水準.05（両側）
x3 <- qPearson(p=.005, N=N, lower.tail=F)	#有意水準.01（両側）
data.frame(x1,x2,x3, row.names=N)		#data.frame()で変数をまとめる


##(b)データの読み込み(男子20人の逸脱行動得点の2時点間の相関係数の検定)
setwd("C:/data/haebara/chap01")
data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T)	


##(c)cor.test(モデル, data, subset)で相関係数の検定
#----------------------------------------------------
#モデル→検定をしたい2変数名を~の後に+でつないで付ける
#data  →検定したいデータが入っているデータセット名を付ける
#subset→男性だけのデータなど，ある一定の基準で被験者を抽出したい場合
#	 変数名 == "質的変数"を入れる
#----------------------------------------------------
cor.test(~中2+小6, data=data1, subset=性別=="男")


##(d)cor.test(x,y)で相関係数の検定
逸脱男 <- subset(data1,性別=="男",select=c(小6,中2))	#男性だけのデータを抽出
attach(逸脱男)						#変数名で参照できるようにする
cor.test(中2,小6)					#2行多く書けば上記と同じことができる


##(e)付表4 t分布の上側確率に対応する値 を再現 p.379
qt(p=0.025,df=18,lower.tail=F)	#qt(p=有意確率,df=自由度,lower.tail=F 上側確率)
qt(p=0.005,df=18,lower.tail=F)	#qt(p=有意確率,df=自由度,lower.tail=F 上側確率)
degree <- c(
	seq(1,20,by=1), seq(22,40,by=2), seq(50,100,by=10), seq(200,400,by=100)
	)			#自由度を指定

x1 <- qt(p=0.05, df=degree,lower.tail=F)	#有意水準.05 （上側）
x2 <- qt(p=0.025,df=degree,lower.tail=F)	#有意水準.025（上側）
x3 <- qt(p=0.005,df=degree,lower.tail=F)	#有意水準.005（上側）
data.frame(x1,x2,x3, row.names=degree)		#data.frame()で変数をまとめる


##(f)t分布を用いた検定(5.8)式を定義

rTOt <- function(r,N){
		r/sqrt(1-r^2) * sqrt(N-2)
			}
rTOt(r=.6,N=20)



################################################
#Step5：検定力の求め方 p.146-47
################################################

##(a) ρ=.4, N=20 をZ変換
atanh(.4)			#4.22式（atanhは逆双曲線正接の関数）
1/sqrt(20-3)			#4.23式


##(b) 棄却の限界値 r=|.444|をz変換
r1 <- qPearson(p=.025, N=20, lower.tail=T)
r2 <- qPearson(p=.025, N=20, lower.tail=F)	#有意水準.05（両側），N=20の場合の
						#棄却の限界値
						
atanh(r1)					#棄却の限界値をz変換
atanh(r2)					#棄却の限界値をz変換


##(c)上下の限界値を平均μzと標準偏差σzで標準化
ZL <- (atanh(r1)-atanh(.4))/(1/sqrt(20-3))
ZU <- (atanh(r2)-atanh(.4))/(1/sqrt(20-3))


##(d)標準正規分布の確率を計算
pnorm(ZL,lower.tail=T)		#下側確率 prob(z< ZL)
pnorm(ZU,lower.tail=F)		#上側確率 prob(z> ZU)
pnorm(ZL,lower.tail=T) + pnorm(ZU,lower.tail=F)	#検出力



################################################
#Step6：相関係数の検定の検定力 （付図1）を再現 p.387
################################################

##(a) 検出力関数を定義 ro=母相関係数, n=サンプルサイズ，sig.level=有意水準
power.cor <- function(ro,n,sig.level){
library(SuppDists)			#SuppDistsライブラリを呼び出す
r1 <-  qPearson(p=sig.level/2, N=n, lower.tail=T)
r2 <-  qPearson(p=sig.level/2, N=n, lower.tail=F)	#棄却域を計算
mu <- atanh(ro)
sigma <- 1/sqrt(n-3)				#roは平均mu, 標準偏差sigmaの正規分布に従う

pnorm(atanh(r1),mean=mu, sd=sigma,lower.tail=T) +
pnorm(atanh(r2),mean=mu, sd=sigma,lower.tail=F)	#検出力
}

power.cor(ro=.4, n=20, sig.level=.05)	#テスト


#(b) 付図1を再現
ro <- seq(-1,1,by=.01)					#母相関係数を作成
power.cor(ro=ro, n=10, sig.level=.05)			#N=10の場合の検定力
plot(power.cor(ro=ro, n=10, sig.level=.05)~ro, type="l",
	xlab="母集団相関係数",ylab="検定力")
lines(power.cor(ro=ro, n=40, sig.level=.05)~ro)		#N=40の場合
lines(power.cor(ro=ro, n=200, sig.level=.05)~ro)	#N=200の場合


################################################
#Step7：区間推定 p.149-
################################################

##(a)任意の相関値を帰無仮説とした検定
qPearson(p=.025, N=20, rho=.4, lower.tail=T)	#H0: ρ=.4の場合の棄却域
qPearson(p=.025, N=20, rho=.4, lower.tail=F)


##(b)図5-5を再現
ro <- seq(-.999,.999,by=.01)		#母集団相関係数を定義
				#(qPearsonはrho=1やrho=-1の場合は使えないので近似値を使う)
qPearson(p=.025, N=20, rho=ro, lower.tail=T)
qPearson(p=.025, N=20, rho=ro, lower.tail=F)

plot(ro~qPearson(p=.025, N=20, rho=ro, lower.tail=T),type="l",
	xlab="標本相関係数",ylab="母集団相関係数")
lines(ro~qPearson(p=.025, N=20, rho=ro, lower.tail=F))


##(c)信頼区間の算出 5.13/14 式
tanh(atanh(.6) - 1.96/sqrt(20-3))
tanh(atanh(.6) + 1.96/sqrt(20-3))
cor.test(~中2+小6, data=data1, subset=性別=="男")	
				#cor.test()でも95percent confidence intervalは算出される
				#cor.test()での信頼区間と5.13/14 式による結果が僅かに違うが
				#元々の標本相関係数が0.5983922と.6と異なるためである



################################################
#Step8：信頼区間を利用したサンプルサイズの決定 p.154-55 & 南風原(1986)
################################################

##(a) 関数の定義（関数の中身は覚えなくても大丈夫）
cor.sample <- function(r,w,conf.level,interval=c(4,10e10)){
	temp <- function(N){
		(tanh(atanh(r)+qnorm((1-conf.level)/2, lower.tail=F)/sqrt(N-3))-
		tanh(atanh(r)-qnorm((1-conf.level)/2, lower.tail=F)/sqrt(N-3)))-w
		}
	ceiling(uniroot(temp,interval)[[1]])
	}

##(b)cor.sample(r=予想される相関係数の値，w=望まれる信頼区間の幅, conf.level=信頼係数)
cor.sample(r=.3, w=.2,conf.level=.95)	#相関係数=.3, 信頼区間の幅=0.2, 95%信頼区間
cor.sample(r=.0, w=.2,conf.level=.95)	#相関係数=0, 信頼区間の幅=0.2, 95%信頼区間

cor.sample(r=.3, w=.2,conf.level=.99)	#相関係数=.3, 信頼区間の幅=0.2, 99%信頼区間
cor.sample(r=.0, w=.2,conf.level=.99)	#相関係数=.0, 信頼区間の幅=0.2, 99%信頼区間