Rで2x3表を網羅的に作ってみる
- 2x3表 *1があるとする。
- これと周辺度数を同じくする表を全部作りたい。。。その数も知りたい。。。
Rでやるなら
makeTables<-function(t){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
x<-seq(0,min(m1[1],m2[1]))
y<-seq(0,min(m1[1],m2[2]))
d11d12<-expand.grid(x=x,y=y)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(data2)
}
- 使い方
t<-c(10,20,30,40,50,60) tables<-makeTables(t)
- 同じく、テーブルの数だけを知りたいとき
countTables<-function(t){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
x<-seq(0,min(m1[1],m2[1]))
y<-seq(0,min(m1[1],m2[2]))
d11d12<-expand.grid(x=x,y=y)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(length(data2[,1]))
}
-使い方
>||
t<-c(10,20,30,40,50,60)
Notables<-countTables(t)
- 周辺度数を算出する関数も使ってます
#行の周辺度数
funcm1<-function(t2){
m1<-c(t2[1,1]+t2[1,2]+t2[1,3],t2[2,1]+t2[2,2]+t2[2,3])
return(m1)
}
#列の周辺度数
funcm2<-function(t2){
m2<-c(t2[1,1]+t2[2,1],t2[1,2]+t2[2,2],t2[1,3]+t2[2,3])
return(m2)
}
makeTablesSkip<-function(t,skip){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
x<-seq(0,min(m1[1],m2[1]),by=skip)
y<-seq(0,min(m1[1],m2[2]),by=skip)
d11d12<-expand.grid(x=x,y=y)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(data2)
}
countTablesSkip<-function(t,skip){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
x<-seq(0,min(m1[1],m2[1]),by=skip)
y<-seq(0,min(m1[1],m2[2]))
d11d12<-expand.grid(x=x,y=y,by=skip)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(length(data2[,1])*skip*skip)
}
- これだと、*2という大きな周辺度数を持つテーブルでも、makeTablesSkip関数では、22281個のテーブルが作られるだけです。そして、テーブルの総数は、おおよそ、2760502500個くらいらしいこともわかります。テーブル数の概数は、桁がずれるほど違うことはないでしょう。
- じゃあ、このskip数を決めるのが面倒くさかったらどうしたらよいでしょう。大まかに言えば、関数の中で動かしている２つの数の掛け算が、「網羅数」に相当しますから、その「網羅しそうな数」に上限を加えてやることで、デフォルトのスキップ数を、入力テーブルから計算させてやると、その目的が果たせそうです。今、10000テーブルを網羅するのが、この我慢の上限のデフォルト値として、場合によっては、もっと我慢強かったり、我慢が足りなかったりするので、それも引数化してやると、次のようになります。
makeTablesSkipDefault<-function(t,limit=10000){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
skip<-ceiling(sqrt((min(m1[1],m2[1])+1)*(min(m1[1],m2[2])+1)/limit))
x<-seq(0,min(m1[1],m2[1]),by=skip)
y<-seq(0,min(m1[1],m2[2]),by=skip)
d11d12<-expand.grid(x=x,y=y)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(list(tables=data2,skip=skip,limit=limit))
}
countTablesSkipDefault<-function(t,limit=10000){
t2<-matrix(t,nrow=2)
m1<-funcm1(t2)
m2<-funcm2(t2)
m<-sum(t2)
skip<-ceiling(sqrt((min(m1[1],m2[1])+1)*(min(m1[1],m2[2])+1)/limit))
x<-seq(0,min(m1[1],m2[1]),by=skip)
y<-seq(0,min(m1[1],m2[2]))
d11d12<-expand.grid(x=x,y=y,by=skip)
d13<-m1[1]-(d11d12[1]+d11d12[2])
d21<-m2[1]-d11d12[1]
d22<-m2[2]-d11d12[2]
d23<-m2[3]-d13
data<-cbind(d11d12[1],d21,d11d12[2],d22,d13,d23)
minval<-apply(data,1,min)
data2<-data[minval>=0,]
return(list(NoTables=length(data2[,1])*skip*skip,skip=skip,limit=limit))
}
- ついでに、HWEな３ディプロタイプから、ずらして2x3表を作ってもみる
sim2x3<-function(af,N1,N2,sd=1){
g<-c(af^2,2*af*(1-af),(1-af)^2)
b<-c(ceiling(g[1]*N1+rnorm(1,mean=N1,sd=sd)),ceiling(g[2]*N1+rnorm(1,mean=N1,sd=sd)),ceiling(g[3]*N1+rnorm(1,mean=N1,sd=sd)),ceiling(g[1]*N2+rnorm(1,mean=N2,sd=sd)),ceiling(g[2]*N2+rnorm(1,mean=N2,sd=sd)),ceiling(g[3]*N2+rnorm(1,mean=N2,sd=sd)))
return(b)
}
#Ntは作成するテーブル数
#N1,N2はケース,コントロール人数
SimTables<-function(N1,N2,sd=1,Nt){
g11<-rep(0,Nt)
g12<-rep(0,Nt)
g13<-rep(0,Nt)
g21<-rep(0,Nt)
g22<-rep(0,Nt)
g23<-rep(0,Nt)
for(i in 1:Nt){
tmp<-sim2x3(runif(1),N1,N2,sd=sd)
g11[i]=tmp[1]
g12[i]=tmp[2]
g13[i]=tmp[3]
g21[i]=tmp[4]
g22[i]=tmp[5]
g23[i]=tmp[6]
}
ret=cbind(g11,g21,g12,g22,g13,g23)
return(ret)
}
