# マルチプルテスティング：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)