2013-07-16

単純な系

kNN 密度推定選択決断

量的計測値がある
その後、量的計測値が対応してとれる
たとえば、治療前の検査値(x)と生存期間(y)、とか
昨日はkNN密度推定を使って、x,yが多次元データであって、しかも、分布が「なんでもあり」なことをやった
今日は思い切り単純にして、xも１次元、yも１次元
ただし、二つの群があることにしよう。この世のことで言えば、「治療Aと治療B」
AとBとには若干の差があるらしいが
個体差も大きく
生存期間の期待値はxの値によらないものの
xの値によって、生存期間のばらつきは大きく影響しているものとする

# knn密度推定関数
knn.density <- function(z,X,k=10){
	d <- length(z)
	tmp <- t(t(X)-z)
	D <- sqrt(apply(tmp^2,1,sum))
	D.s <- sort(D)
	(1:k)/length(X[,1])/(pi^(d/2)/gamma(d/2+1)*D.s[1:k]^d)
}

データを作る

# サンプル数は多い
n <- 1000
# 治療A,Bがx,x2に対応する。治療前の値の分布は同じ
x <- rnorm(n)+2
x2 <- rnorm(n)+2
x <- sort(x)
x2 <- sort(x2)
# 生存期間はA,Bに対してy,y2とする
# x,x2に応じてばらつきが増える。Aの方がBより長い
y <- rep(0,n)
y2 <- rep(0,n)
for(i in 1:n){
	y[i] <- rnorm(1,0,abs(x[i])^0.5*0.5)+2
	y2[i] <- rnorm(1,0,abs(x[i])^0.8*0.6)+1.5
}
plot(c(x,x2),c(y,y2),col=rep(1:2,each=n),pch=20)

同時分布とkNNを使って、xの値(x.check)に対して、yの値(y.check)となる確率・尤度の相対値を計算してみよう

x <- matrix(x,ncol=1)
y <- matrix(y,ncol=1)
x2 <- matrix(x2,ncol=1)
y2 <- matrix(y2,ncol=1)
# 同時分布を作る
XY <- cbind(x,y)
XY2 <- cbind(x2,y2)
# 条件付き確率のグリッド点
x.check <- seq(from=min(c(x,x2)),to=max(c(x,x2)),length=50)
y.check <- seq(from=min(c(y,y2)),to=max(c(y,y2)),length=50)
# 治療 A, Bごとにグリッド点の尤度・確率の相対値を納める
cond.prob <- rep(0,length(y.check))
cond.prob2 <- rep(0,length(y.check))
# kNN法のk
k <- 10
# 求めた相対的尤度・確率から乱数を発生させてそれっぽい分布を得よう
n.r <- 1000
# 乱点を納める行列
rs <- matrix(0,length(x.check),n.r)
rs2 <- matrix(0,length(x.check),n.r)
q.pt0 <- seq(from=0,to=1,by=0.1)
qs <- qs2 <- matrix(0,length(x.check),length(q.pt0))

for(ii in 1:length(x.check)){
	for(i in 1:length(y.check)){
		tmp <- c(x.check[ii],y.check[i])
		cond.prob[i] <- knn.density(tmp,XY,k)[k]
	}
	tmp.p <- cumsum(cond.prob)
	tmp.p <- tmp.p/tmp.p[length(tmp.p)]
	for(j in 1:length(q.pt0)){
		qs[ii,j] <- y.check[which(tmp.p>q.pt0[j])[1]]
	}

	rs[ii,] <- sample(y.check,n.r,replace=TRUE,prob=cond.prob)
	#hist(rs)
	#plot(y,pch=20,col=gray((max(cond.prob)-cond.prob)/(max(cond.prob)-min(cond.prob))))
}
for(ii in 1:length(x.check)){
	for(i in 1:length(y.check)){
		tmp <- c(x.check[ii],y.check[i])
		cond.prob2[i] <- knn.density(tmp,XY2,k)[k]
	}
	tmp.p <- cumsum(cond.prob2)
	tmp.p <- tmp.p/tmp.p[length(tmp.p)]
	for(j in 1:length(q.pt0)){
		qs2[ii,j] <- y.check[which(tmp.p>q.pt0[j])[1]]
	}

	rs2[ii,] <- sample(y.check,n.r,replace=TRUE,prob=cond.prob2)
	#hist(rs)
	#plot(y,pch=20,col=gray((max(cond.prob)-cond.prob)/(max(cond.prob)-min(cond.prob))))
}

治療Aの治療前マーカー値ごとの生存期間をボックスプロットしてみよう。元のデータのばらつきの小さいあたりのボックスは狭く、それ以外は広い

boxplot(t(rs))
boxplot(t(rs2))

２つの治療を比較してみよう。xの値ごとにクオンタイル点をとって、比べる。実線のA、破線のB。有意差はつかないけれど、Aをとる自信の強さには値が付けられそう

q.pt <- c(0.1,0.25,0.5,0.75,0.9)
qqq <- apply(rs,1,quantile,q.pt)
qqq2 <- apply(rs2,1,quantile,q.pt)
matplot(t(rbind(qqq,qqq2)),type="l",col=rep(c(1,2,3,2,1),2),lty=rep(1:2,each=length(q.pt)))
matplot(cbind(qs,qs2),type="l",col=rep(1:10,2),lty=rep(1:2,each=length(q.pt0)))