2013-01-23

２ｘ２表：ぱっと見て感じる

2x2分割表のクイズ
実行すると、２ｘ２表に対応する棒グラフがいくつか並べて提示される
図を先送りすると、そのフィッシャーの正確確率検定ｐ値が表示され、最小ｐ値の棒グラフがハイライトされることで、答え合わせができる
Ｑ『次の３つの棒グラフに対応する２ｘ２表のうち、どれが最小ｐ値をもたらすか？』

# 結構ばらばらな２ｘ２分割表を作る
library(MCMCpack)
# 作成する分割表の数
n.test <- 100
rs <- rdirichlet(n.test,rep(0.9,4))
rs <- round(rs * 50)
# フィッシャーの正確確率検定ですべての表のp値を出す
p <- rep(NA,n.test)
for(i in 1:n.test){
	p[i] <- fisher.test(matrix(rs[i,],2,2))$p
}

# 一度に提示・比較する表の個数を指定
k <- 3
# 『どの表が一番小さいp値をもたらすか？』クイズの実行回数
n.trial <- 50

par(mfcol=c(1,k))
for(i in 1:n.trial){
	# 作成済みの表からk個を選ぶ
	tmp <- sample(1:n.test,k)
	# 表の４セルを２群の赤・黒の棒グラフで示す
	col1 <- c(1,2,0,1,2)
	ylim <- c(0,max(rs[tmp,]))
	for(j in tmp){
		# 『p値は？』
		ttl <- "p=?"
		tmp2 <- rs[j,]
		barplot(c(tmp2[1:2],0,tmp2[3:4]),ylim=ylim,col=col1,main = ttl)
		par(ask=FALSE)
	}
	par(ask=TRUE)
	for(j in tmp){
		# ｐ値を示しつつ
		ttl <- round(p[j],2)
		tmp2 <- rs[j,]
		# 最小ｐ値をもたらした表だけ目立たせる
		if(p[j] == min(p[tmp])){
			barplot(c(tmp2[1:2],0,tmp2[3:4]),ylim=ylim,col=col1,main = ttl)
		}else{
			barplot(c(tmp2[1:2],0,tmp2[3:4]),ylim=ylim,density = rep(100,5),col=col1,main = ttl)
		}
		
		par(ask=FALSE)
	}
	par(ask=TRUE)
}

答え：中央

２ｘ２分割表を別の視覚表現でやってみる

n.iter <- 200
X <- list()
ps <- rep(0,n.iter)
for(i in 1:n.iter){
  n1 <- 300
  n2 <- 200
  p1 <- 0.3
  p2 <- 0.3
  s1 <- sample(0:1,n1,replace=TRUE,prob=c(1-p1,p1))
  s2 <- sample(0:1,n2,replace=TRUE,prob=c(1-p2,p2))
  t <-table(c(s1,s2),c(rep(1,n1),rep(2,n2)))
  X[[i]] <- t
  ps[i] <- chisq.test(t)$p.value  
}
par(ask=TRUE)
for(i in 1:50){
  print(X[[i]])
  ttl <- paste("",X[[i]][1,1],"\t",X[[i]][1,2],"\n",X[[i]][2,1],"\t",X[[i]][2,2],"\n")
  plot(X[[i]],main=ttl)
  #text(c(X[[i]]))
  plot(X[[i]],main=paste("p=",ps[i]))
}
par(ask=FALSE)
hist(ps)
summary(ps)

２群のt検定
- どちらのp値が小さい？

n.iter <- 50*2
X <- list()
ps <- rep(0,n.iter)
for(i in 1:n.iter){
  n1 <- 10
  n2 <- 10
  m1 <- runif(1)*2
  m2 <- runif(1)*2
  s1 <- runif(1)*2
  s2 <- s1
  x1 <- rnorm(n1,m1,s1)
  x2 <- rnorm(n2,m2,s2)
  out <- t.test(x1,x2)
  X[[i]] <- cbind(x1,x2)
  #boxplot(cbind(x1,x2),main=paste("p=",out$p.value))
  ps[i] <- out$p.value
}
par(ask=TRUE)
par(mfcol=c(1,2))
for(i in 1:10){
  par(ask=TRUE)
  boxplot(X[[2*(i-1)+1]],main="p=?")
  boxplot(X[[2*(i-1)+2]],main="p=?")
  par(ask=TRUE)
  boxplot(X[[2*(i-1)+1]],main=paste("p=",ps[2*(i-1)+1]))
  par(ask=FALSE)
  boxplot(X[[2*(i-1)+2]],main=paste("p=",ps[2*(i-1)+2]))
}
par(ask=FALSE)
par(mfcol=c(1,1))
hist(ps)
summary(ps)

- p値を当てよう

n.iter <- 1000
X <- list()
ps <- rep(0,n.iter)
for(i in 1:n.iter){
  n1 <- 10
  n2 <- 10
  m1 <- 0
  m2 <- 2
  s1 <- 2
  s2 <- 2
  x1 <- rnorm(n1,m1,s1)
  x2 <- rnorm(n2,m2,s2)
  out <- t.test(x1,x2)
  X[[i]] <- cbind(x1,x2)
  #boxplot(cbind(x1,x2),main=paste("p=",out$p.value))
  ps[i] <- out$p.value
}
par(ask=TRUE)
for(i in 1:10){
  boxplot(X[[i]],main="p=?")
  boxplot(X[[i]],main=paste("p=",ps[i]))
}
par(ask=FALSE)
hist(ps)
summary(ps)

線形回帰

n.iter <- 100
X <- list()
ps <- rep(0,n.iter)
for(i in 1:n.iter){
  n <- 50
  x <- rnorm(n)
  a <- runif(1)
  b <- runif(1)
  y <- a*x+b
  z <- rnorm(n,runif(1))
  y. <- y + z
  X[[i]] <- cbind(x,y.)
  lm.out <- lm(y.~x)
  s.out <- summary(lm.out)
  ps[i] <- s.out[[4]][2,4]
}
par(ask=TRUE)
for(i in 1:50){
  xlim <- ylim <- range(X[[i]])
  #ttl <- paste("",X[[i]][1,1],"\t",X[[i]][1,2],"\n",X[[i]][2,1],"\t",X[[i]][2,2],"\n")
  plot(X[[i]],xlim=xlim,ylim=ylim)
  #text(c(X[[i]]))
  plot(X[[i]],xlim=xlim,ylim=ylim,main=paste("p=",ps[i]))
}

par(ask=FALSE)
hist(ps)
summary(ps)

次の記事の問題は細胞大小不同