数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

SpherePower()

R Sweave パワー 多次元球 分割表 検定
  • Sweave()を使って解説文書のメモを作る
  • Sweaveのファイル"FxPowerSphere.Rnw"(末尾)を作って
Sweave("FxPowerSphere.Rnw")
  • と実行するとできる"FxPowerSphere.tex"をてふ処理してPDFにしたのがこちら
  • "SpherePower.R"関数はこちら
  • "FxPowerSphere.Rnw"は以下
\documentclass{jarticle}
\author{
Ryo Yamada
}

\title{
関数SpherePowerを説明する文書}

\usepackage[dvips]{graphicx}
\usepackage{graphicx}
\usepackage{bigdelim,multirow}
\usepackage{amsmath,amsthm,amssymb,cases}
\usepackage{ascmac}
\usepackage{url}
\usepackage{eclbkbox}
\usepackage{wrapfig}
\usepackage{listings, jlisting}
\usepackage[dvips,usenames]{color}
 \usepackage{makeidx}
\usepackage{url}
\usepackage[dvipdfm,%
bookmarks=true,bookmarksnumbered=true,bookmarkstype=toc,%
colorlinks=true,% 
%colorlinks=false,%
linkcolor=blue,%
citecolor=blue,%
filecolor=blue,%
pagecolor=blue,%
urlcolor=blue%
]{hyperref}
\renewcommand{\baselinestretch}{1}
\setlength{\textheight}{33\baselineskip}


\makeindex


\begin{document}
\maketitle
\lstset{%
 language={R},
 backgroundcolor={\color[gray]{.85}},%
 basicstyle={\small\tt },%
% basicstyle={\small },%
 identifierstyle={\small},%
 commentstyle={\small\itshape},%
 keywordstyle={\small\bfseries},%
 ndkeywordstyle={\small},%
 stringstyle={\small\ttfamily},
 frame={tb},
 breaklines=true,
 columns=[l]{fullflexible},%
 numbers=left,%
 xrightmargin=0zw,%
 xleftmargin=0zw,%
% xleftmargin=3zw,%
 numberstyle={\scriptsize},%
 stepnumber=1,
 numbersep=1zw,%
 lineskip=-0.5ex%
}
\section{"SpherePower()"関数}
\begin{lstlisting}
ソースをまず示す
SpherePower<-function(df,Kv,Ks,CheckDists,Tiis,nperm=1000){
	Ks[which(Ks==0)]<-0.0001
	Kss<-matrix(rep(Kv,length(Ks)),nrow=length(Ks),byrow=TRUE)*Ks
	Ntest<-length(Tiis[,1])
	rss<-RandomSphere(df=df,n=nperm)
	PowerOut<-matrix(0,length(Ks),length(CheckDists))
	for(ik in 1:length(Ks)){
		Pii<-Kss[ik,]
		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 <- NULL
			if (length(SerialStat) > 1) {
				SStat1 <- SerialStat[1:(length(SerialStat) - 1)]
				SStat2 <- SerialStat[2:length(SerialStat)]
				intMyX <- outer(SStat1, CheckDists, FUN = "-") * outer(SStat2, CheckDists, FUN = "-") # 線分が、CheckDistsの値を横切っていれば(-1)、同じ側なら(+1)
			}
			
			intMyX <- rbind(intMyX, (SerialStat[length(SerialStat)] - CheckDists) * As[SerialTopT[length(SerialTopT)]])
			crossID<-which(intMyX<0,arr.ind=TRUE)
			noncrossID<-which(intMyX>=0,arr.ind=TRUE)
			intMyX[crossID]<--1
			intMyX[noncrossID]<-0
			zeros <- which(apply(abs(intMyX), 2, sum) == 0) # 交叉しないCheckDistsID？# 20110327
			intMyX3 <- -intMyX * sign(As[SerialTopT]) # 統計量の籠の外側から内側へ入るIDは負、内側から外側へ出るIDは正
			intMyX2 <- intMyX^apply((-intMyX), 2, cumsum) * intMyX
			# 交点を出す
			Crossings <- t(outer(CheckDists, Bs[SerialTopT], FUN = "-"))/As[SerialTopT]
			ylim <- c(min(Crossings), max(Crossings))
			ChiProbs <- pchisq(Crossings^2, df, lower.tail = FALSE)
			AnswerPs <- apply(intMyX3 * ChiProbs, 2, sum)
			# Inside-Outside, Outside-Inside-Outsideのパターン別の計算
			# 0の扱い AnswerPsが実質的に0のとき、1を与えるか0を与えるかは
			# スタート地点が内側か外側かで決まる
			AnswerPs[which(AnswerPs < 0)] <- 1 + AnswerPs[which(AnswerPs < 0)]
			AnswerPs[which(AnswerPs > 0)] <- AnswerPs[which(AnswerPs > 0)]
			outside<-which((SerialStat[1]-CheckDists)>0 & AnswerPs==0)
			AnswerPs[outside] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 1)/2
			AnswerPs[zeros] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 1)/2
			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)
}
\end{lstlisting}
順次、解説して行こう
\subsection{引数について}
\begin{itemize}
\item $df$:自由度・次元。分割表検定なので、自由度がある。これは、標準化正単体座標系では次元である。
\item $Kv,Ks$:対立仮説は、自由度次元の空間の点で定義する。
この対立仮説の点を、その方向の単位ベクトル$Kv$と長さ$Ks$で指定する。
$Ks$は0以上の実数ベクトルとして、
同一方向の複数の仮説を一括して指定することができる。
$Ks^2$は、自由度$df$のピアソンの独立性検定を行ったときの、
$\chi^2$の値の最尤値である。
今、$N\times M$表が与えられているときに、それに対応する$Kv,Ks$を"RegularSpherize()"関数で算出
し、それを"SpherePower()"関数の引数として与えることもできる
\item $Tiis$:自由度１のテストに対応する単位法線ベクトルを１つ以上実施するとき、単位法線ベクトルを
行ベクトルとする、行列が$Tiis$である。
テストの単位法線ベクトルは、向きのあるテストを表すから、両側テストの場合には、単位ベクトルと、
その反対向きの単位ベクトルの２ベクトルをT$Tiis$に含めること。
この$Tiis$は座標空間で等高線の形を定める。
等高線の形は、自由度次元球に接する平面が作る凸包である。
$N\times M$表のセルに重みづけを与えることで自由度１のテストを定義する場合には
"TestDf1()"を用いて、その自由度次元ベクトルを算出することができる。
\item $CheckDists$:帰無仮説検定において、ある$\chi^2$値を検定統計量とするとき、その値の平方根を
$CheckDists$として与える。
複数の値をリストとして与えることができる。
この値は、座標空間で、等高線の位置を定める。
個々の自由度１テストで$CheckDists^2$の$\chi^2$の値に相当する平面がこの等高線であり、
半径$CheckDist$の$df$次元球に接する凸包である。
\item $nperm$: 検出力を算出するにあたって、発生させる乱方向の数
\end{itemize}
\subsection{検定統計量の等高線}
自由度２の場合に、自由度１の両側検定を$Nt$個実施するとし、検定統計量$Kt^2$の等高線を描く。
自由度２の$\chi^2=Kt^2$の等高線は円を成すが、それを黒で、自由度１の検定量の直線を赤で描く。
複数の自由度１のテストのうちの最大の統計量を採用する検定では、
この直線のもっとも内側の線を結んだ凸包が等高線となる。
<<fig=TRUE>>=
library(sphere)
df<-2
Nt<-4
rt<-runif(Nt)*2*pi
Tiis<-cbind(cos(rt),sin(rt))
Tiis<-rbind(Tiis,-Tiis)
t<-seq(from=0,to=1,length.out=100)*2*pi
Kt<-runif(1)*5
plot(cos(t)*Kt,sin(t)*Kt,xlim=c(-Kt*2,Kt*2),ylim=c(-Kt*2,Kt*2),type="l")
for(i in 1:length(Tiis[,1])){
	abline(1/Tiis[i,2]*Kt,-Tiis[i,1]/Tiis[i,2],col=2)
}
@
\subsection{帰無仮説が真であるときの自由度次元標準正規分布}
帰無仮説が真であるとき、
そのような集団から標本をサンプリングするとき、
得られる標本はその点を中心とした、自由度次元標準正規分布を取る。
そのような分布をプロットしてみる。
前の等高線に重ねて描いてみる。
<<fig=TRUE>>=
plot(cos(t)*Kt,sin(t)*Kt,xlim=c(-Kt*2,Kt*2),ylim=c(-Kt*2,Kt*2),type="l")
for(i in 1:length(Tiis[,1])){
	abline(1/Tiis[i,2]*Kt,-Tiis[i,1]/Tiis[i,2],col=2)
}
Ns<-10000
X<-matrix(rnorm(Ns*df),ncol=df)
par(new=TRUE)
plot(X[,1],X[,2],cex=0.01,xlim=c(-Kt*2,Kt*2),ylim=c(-Kt*2,Kt*2),col=3)
@
\subsection{対立仮説が真であるときの自由度次元標準正規分布}
対立仮説が真であるとき、その点は$Ks\times Kv$。
そのような集団から標本をサンプリングするとき、
得られる標本はその点を中心とした、自由度次元標準正規分布を取る。
そのような分布をプロットしてみる。
前の等高線に重ねて描いてみる。
<<fig=TRUE>>=
Ns<-10000
X<-matrix(rnorm(Ns*df),ncol=df)
Ks<-runif(1)*5
Kv<-rnorm(df)
Kv<-Kv/sqrt(sum(Kv^2))
X[,1]<-X[,1]+Ks*Kv[1]
X[,2]<-X[,2]+Ks*Kv[2]
a<-max(abs(X))
xlim<-c(-a,a)
plot(cos(t)*Kt,sin(t)*Kt,xlim=xlim,ylim=xlim,type="l")
for(i in 1:length(Tiis[,1])){
	abline(1/Tiis[i,2]*Kt,-Tiis[i,1]/Tiis[i,2],col=2)
}

par(new=TRUE)
plot(X[,1],X[,2],cex=0.01,xlim=xlim,ylim=xlim,col=3)
@
\subsection{真値からある方向へずれる場合}
\subsubsection{あるテストの統計量}
今、ある点$Ks\times Kv$から、単位ベクトル$x$の方向へとずれる場合を
考える。

このずれた点を$Ks\times Kv + u x$とする。

また、ある自由度１のテストの法線単位ベクトルが$t_i$であるとき、
その統計量は
内積
$$<Ks\times Kv + u x, t_i>=Ks<Kv,t_i> + u <v,t_i>$$
で算出できる。
これは、統計量が、始点からの距離$u$の一次関数であることを示している。
以下の図は、ある点からの一定方向へのずれと、あるテストの等高線とを描いたものである。

<<fig=TRUE>>=
plot(cos(t)*Kt,sin(t)*Kt,xlim=xlim,ylim=xlim,type="l")
Kts<-seq(from=0,to=1,length.out=100)*Kt*3
i<-1
for(j in 1:length(Kts)){
	abline(1/Tiis[i,2]*Kts[j],-Tiis[i,1]/Tiis[i,2],col=gray(j/length(Kts)))
	abline(-1/Tiis[i,2]*Kts[j],-Tiis[i,1]/Tiis[i,2],col=gray(j/length(Kts)))
}
par(new=TRUE)
plot(Ks*Kv[1],Ks*Kv[2],pch=19,cex=2,col=3,xlim=xlim,ylim=xlim)
rs<-RandomSphere(df,1)
u<-seq(from=0,to=10,length.out=100)
par(new=TRUE)
plot(Ks*Kv[1]+rs[1]*u,Ks*Kv[2]+rs[2]*u,type="l",col=3,xlim=xlim,ylim=xlim)
@
\subsubsection{複数のテスト}
今、$Nt$個の自由度１の両側テストを行うとき、それぞれの統計量はある真の点からの距離の
一次関数になる。それを横軸に真の点からの距離、縦軸に統計量で示せば、以下のようになる。
複数のテストを実施し、その統計量の最大値を採用することにすれば、複数の直線が作る上縁を成す
折れ線が統計量である。
<<fig=TRUE>>=
df<-2
t<-runif(1)*2*pi
Kv<-c(cos(t),sin(t))
Ks<-sample(1:50,1)
CheckDists<-seq(from=0,to=sqrt(Ks^2+df-1)*1.2,length.out=50)

Nt<-20
library(MCMCpack)
Tiis<-RandomSphere(df=df,n=Nt)
Tiis<-rbind(Tiis,-Tiis)
nperm<-100
	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))
ik<-1
i<-1
		Pii<-Kss[ik,]
		# nperm個の方向について、指定検定統計量より大なる統計量を得る確率を格納する
		AccumProbs <- matrix(0, nperm, length(CheckDists))
			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)
			xlim<-c(-max(abs(Bs))*2,max(abs(Bs))*2)
plot(0,0,pch=19,cex=1,xlim=xlim,ylim=xlim)
for(u in 1:length(Tiis[,1])){
	abline(Bs[u],As[u],col=2)
}
abline(v=0)
abline(h=0)
@
\section{パワー}
\subsection{検定統計量がある値より大きいとき}
自由度１のテストを複数実施し、
その中で最大の統計量を採用することは、
原点を中心とした球に接する平面が作る凸包状の等高線によって統計量を定めることであることを示してきた。

また、ある真の点から、ある方向へずれていく場合、そのずれの大きさと統計量との関係は、一次関数で表され、
複数のテストの最大統計量は、折れ線グラフになることも示した。

今、統計検定量が与えられると、その値以上の値を取るか否かは、真の点からの距離で定まるが、それは、グラフで表せば、
以下のように、検定統計量を表す直線が水平線で示され、折れ線がそれよりも上部にあれば、統計量は
検定統計量よりも大きく、下部にあれば、検定統計量よりも小さいことを示す。
<<fig=TRUE>>=
plot(0,0,pch=19,cex=1,xlim=xlim,ylim=xlim)
for(u in 1:length(Tiis[,1])){
	abline(Bs[u],As[u],col=3)
}
abline(v=0)
abline(h=0)
abline(h=max(Bs)/2,col=2)
@
\subsection{パワーは検定統計量以上の統計量を観測する確率}
\subsubsection{ある方向へのずれの場合}
対立仮説が真であるときに、検定統計量よりも大きい統計量を観察する確率が
パワーの定義である。
したがって、真の点から、特定の方向にずれる場合のみを考えれば、統計量の
折れ線(緑の線の上縁)が、閾値統計量の線(赤の水平線)の上部に来ているような、ずれを起こす確率
がパワーである。
<<fig=TRUE>>=
df<-2
t<-seq(from=0,to=1,length.out=100)*2*pi
Kt<-runif(1)*5
plot(cos(t)*Kt,sin(t)*Kt,xlim=c(-Kt*2,Kt*2),ylim=c(-Kt*2,Kt*2),type="l")
slope1<-runif(1)
slope2<-slope1+0.1
abline(0,slope1)
abline(0,slope2)
@
今、この円(球)の外側の点を観測する確率は、自由度$df$のカイ二乗分布に従う。
球の内外の確率比はカイ二乗分布で確定するとも言える。
ここで、
この図に示したように、ある楔状の部分を考える。
この部分は、空間全体を原点からの距離について均等に縮めたものであるから、
この楔状の部分における、球の内外の確率比は、球全体のそれと等しい。

\subsubsection{特定方向の統計量の折れ線グラフと検定統計量の水平線の交叉関係を算出するアルゴリズム}
特定の方向において、統計量の折れ線と、検定統計量の水平線との交叉関係を次のように求めることとする。
\begin{itemize}
\item 距離0において、すべてのテストの統計量を算出し、その最大値を与えるテストを確認する
\item 移動方向について、上縁を形成する直線を探索する
\item 上縁を形成する直線が、ある直線から次の直線に交代するときには、交代前後の２直線の交点にて交代する。
したがって、上縁を構成している直線が与えられたとき、移動方向に交点を有する直線が、その候補である
\item その候補直線のうち、もっとも少ない移動量で交叉する直線が、次に上縁を構成する直線である
\item また、ひとたび、候補から漏れた直線は、再度候補になることはない
\item このことを利用して、候補直線がなくなるまで上縁を成す直線を探索することで、上縁が得られる
\item 上縁の折れ線のうち、水平な直線(検定統計量)と交叉する線分は、
線分の両端の点の垂直座標(統計量)のうちの一つは検定統計量より大きく、もう一つが小さいものであることで確認できる
\item 交点を持つ線分を含む直線と水平直線との交点を求めることで、折れ線と水平直線との交点が列挙できる
\item また、その折れ線線分が水平直線の下から上へよぎるか、上から下へよぎるかは、線分の傾きで決まることから、
その情報も得られる
\item 交点の始点からの距離がわかるので、その遠近の確率比は自由度$df$のカイ二乗分布に従うことから、折れ線のうち、
検定統計量の上部をなす範囲の確率が算出できる
\end{itemize}
\subsection{すべての方向は平等なので、乱方向について検討し、その平均を求める}
特定の方向について、検定統計量以上を観察する確率が求められたから、十分に多い方向をランダムに
発生させ、そのそれぞれの方向について、検定統計量以上を観測する確率を求め、
全方向に関して、平均を求めれば、それがパワーである。
\end{document}

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------