マルチプルテスティング：James-Stein推定

• こちらから

my.JS.estimate <- function(x){
	m <- mean(x)
	A1 <- (length(x)-2)/sum((x-m)^2)
	y <- (1-A1)*(x-m) +m
	return(y)
}
n.iter <- 1000
N <- 1000
A <-pi
R.As <- R.Bs <- rep(0,n.iter)
for(i in 1:n.iter){
	# 変なデータとして、正規分布と指数分布の合わさったものをつくってみる
	N <- 100
	N2 <- 50
	A <- 4
	A2 <- 6
	mu <- rnorm(N,3,sqrt(A))
	mu2 <- rexp(N,1)+3
	z.true <- c(mu,mu2)
	z <- mu + rnorm(N,0,1)
	z2 <- mu2 + rnorm(N,0,2)
	z.obs <- c(z,z2)
	z.JS <- my.JS.estimate(z.obs)
	R.As[i] <- sum((z.obs-z.true)^2)
	R.Bs[i] <- sum((z.JS-z.true)^2)
}
xlim <- ylim <- range(c(R.As,R.Bs))
plot(R.As,R.Bs,xlim = xlim,ylim=ylim)
abline(0,1,col=2)

• 真値が0(帰無仮説が真のz検定)でもJS統計量のリスク関数値は小さくなる

n.iter <- 1000
N <- 1000
A <-pi
R.As <- R.Bs <- rep(0,n.iter)
for(i in 1:n.iter){
	N <- 100
	A <- 4
	mu <- rep(0,N)
	z.true <- c(mu)
	z <- mu + rnorm(N,0,1)
	z.obs <- z
	z.JS <- my.JS.estimate(z.obs)
	R.As[i] <- sum((z.obs-z.true)^2)
	R.Bs[i] <- sum((z.JS-z.true)^2)
}
xlim <- ylim <- range(c(R.As,R.Bs))
plot(R.As,R.Bs,xlim = xlim,ylim=ylim)
abline(0,1,col=2)

• カイ二乗検定でも

n.marker <- 100
n.case <- n.cont <- 100
n.iter <- 100
phenotype <- rep(0:1,n.case)
X <- matrix(sample(0:1,n.marker*length(phenotype),replace=TRUE),ncol=n.marker)
R.As <- R.Bs <- rep(0,n.iter)
for(i in 1:n.iter){
	sh.phenotype <- sample(phenotype)
	chis.list <- apply(X,2,chisq.test,sh.phenotype)
	chis <- rep(0,n.marker)
	for(j in 1:n.marker){
		chis[j] <- chis.list[[j]][[1]]
	}
	chis.js <- my.JS.estimate(chis)
	R.As[i] <- sum((chis-0)^2)
	R.Bs[i] <- sum((chis.js-0)^2)

}
xlim <- ylim <- range(c(R.As,R.Bs))
plot(R.As,R.Bs,xlim = xlim,ylim=ylim)
abline(0,1,col=2)

• LDを入れると統計量のバリエーションが減るので、JS推定量はさらに圧縮される

n.marker <- 100
n.case <- n.cont <- 100
n.iter <- 100
phenotype <- rep(0:1,n.case)
X <- matrix(sample(0:1,n.marker*length(phenotype),replace=TRUE),ncol=n.marker)
R.As <- R.Bs <- rep(0,n.iter)
R.As.LD <- R.Bs.LD <- rep(0,n.iter)
for(i in 1:n.iter){
	sh.phenotype <- sample(phenotype)
	chis.list <- apply(X,2,chisq.test,sh.phenotype)
	chis <- rep(0,n.marker)
	for(j in 1:n.marker){
		chis[j] <- chis.list[[j]][[1]]
	}
	chis.js <- my.JS.estimate(chis)
	# 全マーカーを10グループに分けて、１グループ内は「absolute LD」
	chis.LD <- rep(chis[1:10],10)
	chis.LD.js <- my.JS.estimate(chis.LD)

	R.As[i] <- sum((chis-0)^2)
	R.Bs[i] <- sum((chis.js-0)^2)
	R.As.LD[i] <- sum((chis.LD-0)^2)
	R.Bs.LD[i] <- sum((chis.LD.js-0)^2)

}
xlim <- ylim <- range(c(R.As,R.Bs))
plot(R.As,R.Bs,xlim = xlim,ylim=ylim)
points(R.As.LD,R.Bs.LD,col=3)
abline(0,1,col=2)

plot(sort(chis.LD.js))
points(1:n.marker,sort(chis.js),col=2)