モンテカルロで自由度１の多重検定検出力

sphere MoteCarloPowerChisq 多次元球 分割表
  • Sweaveを用いて作成した解説メモはこちら
  • 予定パッケージ名
  • 関数名
    • SpherePower
  • タイトル
SpherePower
  • 説明
真値座標を指定、その指定倍距離の点を真値とする。
検定閾値距離を指定。
自由度１テストの法線ベクトルのセットを指定。
その最大値を統計量とする、多重検定を実施しているシチュエーション

そのうえで、モンテカルロ法にて、検定閾値を超える統計量を観察する。
  • 使用例
SPout<-SpherePower(df=df,Kv=Kv,Ks=Ks,CheckStats=CheckStats,Tiis=Tiis,nperm=nperm)
  • ソース
SpherePower<-function(df,Kv,Ks,CheckDists,Tiis,nperm=1000){
	Ks[which(Ks==0)]<-0.0001 # Ks=0は0から少しずらしておく。計算の都合による。
	# 同一方向、指定長さのベクトルのセットを作る
	Kss<-matrix(rep(Kv,length(Ks)),nrow=length(Ks),byrow=TRUE)*Ks
	# テストベクトルの数
	Ntest<-length(Tiis[,1])
	# 乱方向ベクトルのセットを発生する
	rss<-RandomSphere(df=df,n=nperm)
	# 出力格納ベクトル
	# Ksの値数ｘ検定統計量の値数のパワー値を算出する
	PowerOut<-matrix(0,length(Ks),length(CheckDists))
	# Ksの値ごとに計算を繰り返す
	for(ik in 1:length(Ks)){
		Pii<-Kss[ik,]
		# nperm個の方向について、指定検定統計量より大なる統計量を得る確率を格納する
		AccumProbs <- matrix(0, nperm, length(CheckDists))
		# 乱方向ごとに計算する
		for (i in 1:nperm) {
			rs <- matrix(rss[i,],nrow=1)
			tmptheta <- acos(cosxy(t(Pii), rs))
			C <- nrm(Pii) * cos(tmptheta)
			D <- nrm(Pii) * sin(tmptheta)
			As <- t(Tiis %*% t(rs))
			Bs <- t(Tiis %*% Pii)
			current <- which(Bs == max(Bs))[1]
			SerialTopT <- c(current)
			currentR <- 0
			SerialR <- c(currentR)
			SerialStat <- c(Bs[current])
			CheckList <- 1:Ntest
			CheckList <- CheckList[which(CheckList != current)]
			while (length(CheckList) > 0) {
				tmpks <- -(Bs[current] - Bs[CheckList])/(As[current] - 
				As[CheckList])
				tmpList <- which(tmpks > currentR)
				if (length(tmpList) > 0) {
					current <- CheckList[which(tmpks == min(tmpks[tmpList]))][1] # 次の等高線を作る線分のID
					currentR <- min(tmpks[tmpList]) # 対立仮説点からの距離
					SerialTopT <- c(SerialTopT, current) # IDをリストに
					SerialR <- c(SerialR, currentR) # 対立仮説点からの距離をリストに
					SerialStat <- c(SerialStat, As[current] * currentR + Bs[current]) # 統計量のリスト
				}	
				CheckList <- CheckList[tmpList] # 検討するべきIDのリストを更新
				CheckList <- CheckList[which(CheckList != current)] # 今登録したIDを除く
			}
			# intMyX:交叉すれば1
			# intMyX3:外から内へ入る場合に-1,内から外へ出る場合に+1
			# 線分は半開区間で、その始点は含み、終点は含まない
			intMyX<-intMyX3 <- matrix(0,length(CheckDists),length(SerialTopT))
			AnswerPs<-rep(0,length(CheckDists))
			for(j in 1:length(CheckDists)){
				tmpSerialStat<-SerialStat-CheckDists[j]
				#tmpSerialStat[which(abs(tmpSerialStat)<10^(-10))]<-0
				#print(tmpSerialStat)
				if(tmpSerialStat[1]>=0){
					AnswerPs[j]<-1
				}
#print("pre")
#print(AnswerPs)
				if(length(tmpSerialStat)>1){
					for(jj in 1:(length(tmpSerialStat)-1)){
						if(tmpSerialStat[jj]==0){
							intMyX[j,jj]<-1
							intMyX3[j,jj]<-sign(As[SerialTopT[jj]])
						}else{
							if(tmpSerialStat[jj]*tmpSerialStat[jj+1]<0){
								intMyX[j,jj]<-1
								intMyX3[j,jj]<-sign(As[SerialTopT[jj]])
							}
						}
					}
				}
				
				if(tmpSerialStat[length(tmpSerialStat)]==0){
					intMyX[j,length(tmpSerialStat)]<-1
					intMyX3[j,length(tmpSerialStat)]<-sign(As[SerialTopT[length(tmpSerialStat)]])
				}else{
					if(tmpSerialStat[length(tmpSerialStat)]*As[SerialTopT[length(tmpSerialStat)]]<0){
						intMyX[j,length(tmpSerialStat)]<-1
						intMyX3[j,length(tmpSerialStat)]<-sign(As[SerialTopT[length(tmpSerialStat)]])
					}
				}
			}


			zeros <- which(apply(abs(intMyX), 1, sum) == 0) # 交叉しないCheckDistsID？
			Crossings <- t(outer(CheckDists, Bs[SerialTopT], FUN = "-"))/As[SerialTopT]
			#print(intMyX3)
			#print(intMyX)
			#print(zeros)
			ChiProbs <- pchisq(Crossings^2, df, lower.tail = FALSE)
			#print(Crossings^2)
			#print(ChiProbs)
			AnswerPs <- AnswerPs + apply(t(intMyX3) * ChiProbs, 2, sum)
			AnswerPs[zeros]<-AnswerPs[zeros] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 1)/2
			#print(AnswerPs)
			AnswerPs[which(AnswerPs<0)]<-1+AnswerPs[which(AnswerPs<0)]
			AnswerPs[which(AnswerPs>1)]<-AnswerPs[which(AnswerPs>1)]-1
			AccumProbs[i, ] <- AnswerPs
		}
		PowOut <- apply(AccumProbs, 2, sum)/nperm
		PowerOut[ik,]<-PowOut
	}

	list(df=df,Kv=Kv,Ks=Ks,CheckDists=CheckDists,Tiis=Tiis,nperm=nperm,p.out=PowerOut)
}
  • Rdファイル
\name{SpherePower}
\alias{SpherePower}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
SpherePower
}
\description{
真値座標を指定、その指定倍距離の点を真値とする。
検定閾値距離を指定。
自由度１テストの法線ベクトルのセットを指定。
その最大値を統計量とする、多重検定を実施しているシチュエーション

そのうえで、モンテカルロ法にて、検定閾値を超える統計量を観察する。
}
\usage{
SPout<-SpherePower(df=df,Kv=Kv,Ks=Ks,CheckStats=CheckStats,Tiis=Tiis,nperm=nperm)
}
\arguments{
\item{df}{自由度・次元}
\item{Kv}{真値方向ベクトル、単位ベクトルなら解釈が容易}
\item{Ks}{真値の原点からの距離}
\item{CheckDists}{検定の統計量(原点からの距離)。この二乗がカイ二乗値}
\item{Tiist}{法線ベクトルを行ベクトルとする行列}
\item{nperm}{モンテカルロ乱点数}
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
\item{df}{自由度・次元}
\item{Kv}{真値方向ベクトル、単位ベクトルなら解釈が容易}
\item{Ks}{真値の原点からの距離}
\item{CheckDists}{検定の統計量(原点からの距離)。この二乗がカイ二乗値}
\item{Tiist}{法線ベクトルを行ベクトルとする行列}
\item{nperm}{モンテカルロ乱点数}
\item{p.out}{Ks値を行に、CheckDistsを列として、CheckDists半径球の外側に点を得る確率}
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}
\examples{
df<-sample(1:10,1)
t<-runif(1)*2*pi
#Kv<-c(cos(t),sin(t),0)
Kv<-runif(df)
Kv<-Kv/sqrt(sum(Kv^2))
Ks<-sample(1:50,1)
CheckDists<-seq(from=0,to=sqrt(Ks^2+df-1)*1.2,length.out=50)
#CheckDists<-Ks+df-1

#Tiis<-matrix(c(1,0,0,0,1,0,0,0,1),3,3,byrow=TRUE)
Nt<-1
Nt<-sample(1:1000,1)
library(MCMCpack)
Tiis<-rdirichlet(Nt,rep(1,df))
Tiis<-rbind(Tiis,-Tiis)
nperm<-100
spout<-SpherePower(df=df,Kv=Kv,Ks=Ks,CheckDists=CheckDists,Tiis=Tiis,nperm=nperm)

plot(c(spout$p.out),type="l")
par(new=TRUE)
plot(pchisq(CheckDists^2,df,Ks^2,lower.tail=FALSE),col=2,type="l")