2010-10-05

疾患のモデル化2

ryamada222010-10-05

2010-10-05

どこで区切りたいのか(診断基準をデータ処理的に一般化する)

R 関節リウマチ

(昨日)は、関節リウマチの多関節炎・対称性関節炎という特徴が、相互にどういう関係にあるかということをとっかかりに、ポアッソン分布を少しいじった
少し進めてみるとともに、一般化することとする
ある個人tは関節リウマチという病的全身状態にある
- 右と左に炎症関節をそれぞれ $n_{t,1},n_{t,2}$ 個持つとする
- 今、「左右」には特に意味がないとする
- このとき $n_{t,1}+n_{t,2}$ が炎症関節総数であって、その値は0であることは少なく(関節外症状先行型で関節炎未顕段階など)、複数の関節が炎症を持つことが多い。その数は当然離散的で、数えられる程度の値であり、多すぎない値が現実的なので、この $n_{t,1}+n_{t,2}$ はポアッソン分布に従って、その平均値は、２、３、４など、その位とみなしてよいのではないだろうか。実際には、「平均値」を統計的に知ることももちろんできる
- 多関節炎かどうかの判断は $n_{t,1}+n_{t,2} \ge 2$ と数式ではあらわされる
- それでは、対称性関節炎はどうなるだろうか。 $min(n_{t,1},n_{t,2}) \ge 1$ となる。
- 対称性関節炎条件は、明らかに、多関節炎条件よりも厳しい条件で、それは、対称性関節炎条件は多関節炎条件に含まれるとも言い換えられる。
これをもう少し一般化しよう
- k個の群がある。リウマチの場合には、右と左の２群だったものとする
- それぞれの群は、0,1,2,...,という離散的な量的値をとるものとする。リウマチの例では、これが、右(もしくは)左の炎症性関節数に相当する。
- このようなとき、個人tの状態は $N_t=\{n_{t,1},n_{t,2},...,n_{t,k}\}$ という長さkのベクトルによってあらわされる
- ここで、●かそうでないかという２分診断をすることを考える
- それは、k次元空間を２つの部分空間に分ける面を定義することである
- リウマチの多関節炎定義や対称性関節炎定義を、このk次元空間に持ち込んで、さらに、閾値についても一般化することとする
- 多関節炎定義の一般化
  - $\sum_{i=1}^k n_{t,i} \ge K$ 、 $K$ は閾値。多関節炎の場合は $K=2$
- 対称性関節炎定義の一般化
  - $\forall\;i n_{t,i} \ge K$ (対称性の場合は $K=1$ )
  - 逆もあって $\exists i \;n_{t,i} \ge K$ は容易に考えられるもの
  - さらに、「すべての」「ひとつでも」をもう少し、量的に扱えば、「k個のうち、 $l$ 個の $i$ にあって、 $n_{t,i} \ge K$ であり残りの $k-l$ 個の $i$ では[tex:n_{t,i}
- そのほかには、適当にどうにでも区切りを考えて入れることができる

2010-10-05

ポアッソン分布の分解

R 関節リウマチ

さて、炎症性関節総数がポアッソン分布に従うときに、右と左の炎症性関節数はどんな分布になるのだろうか
それに進む前に、多変量ポアッソン分布についてリンクを張っておく。多変量正規分布があるように、複数の軸について、それぞれポアッソン分布に従うような多変量を扱うことができて、それが多変量ポアッソン分布(Mathematicaのページがこちら)
では、単変量ポアッソン分布に従うような炎症関節数があるときに、左右の振り分けが独立であったとしたら、左右のその数は、どんな風になるだろうか。分布は、ポアッソン分布と似た感じの非対称な離散的１峰性分布になる。それを確かめるのが次のRのソース。ただし、k次元に一般化して処理している。次元数(k)を２にしておけば、左右の話しになる。

k<-2 # 次元 変数の数
n<-30 # 考慮する、次元ごとの値の最大値
npoly<-6 # 全次元について足し合わせ、単変量ポアッソン分布を仮定したときの平均値
ttt<-dpois(0:n,npoly) # 単変量ポアッソン分布の確率密度分布

aaa<-array(0,rep((n+1),k)) # 計算結果格納用アレイ
keta<-0:(k-1) # n進数計算のための道具
keta<-(n+1)^keta # n進数のアレイ内番地計算用の道具

counter<-rep(0,k) # 同上
counter[1]<--1
loop<-TRUE
while(loop){
	counter[1]<-counter[1]+1
	for(i in 1:(k-1)){
		if(counter[i]==n+2){
			counter[i]<-0
			counter[i+1]<-counter[i+1]+1
		}
	}
	banchi<-sum(keta*counter)+1
	vvv<-counter+1
	vvv2<-sum(vvv)-k+1
	if(vvv2<=n+1){
		aaa[banchi]<-ttt[vvv2]*exp(-sum(lgamma(vvv)))
	}
	if(sum(counter)==(n+1)*k){
		loop<-FALSE
	}
}
par(mfcol=c(2,2))
#image(aaa[,,1]) # k=3の場合
image(aaa)
plot(apply(aaa,1,sum),type="b",main="軸１の分布")
plot(apply(aaa,2,sum),type="b",main="軸２の分布")


ax<-apply(aaa,1,sum)
ax<-ax/sum(ax)
tmpmean<-sum((0:n)*ax)
plot(0:n,dpois(0:n,tmpmean),ylim=c(0,max(ax)))

par(new=TRUE)
plot(0:n,ax,type="l",col="red",ylim=c(0,max(ax)),main="1軸の分布に、その平均値をパラメタとして与えたポアッソン分布を重ねた")