ガンマ分布は、期間 ごとに1回くらい起こるランダムな事象が
回起こるまでの時間の分布としてモデル化されることがあるらしい。
ガンマ分布自体は、shape とscale (もしくは逆数のrate) のふたつのパラメータで決定される。
悪性腫瘍(いっぱんにがん)は、遺伝子に変異が蓄積してがんとして発病する、と言われる。ここで、ひとりのヒトに変異が1年あたり 回入るとき、これが
回蓄積すると発癌する、とする。ここで、変異の修復は考えない。
1年あたりに1回以下の珍しい頻度なので、変異が入る事自体はポアソン分布から得られる。
ベタにシミュレーションしつつ、ある年で生じた変異を翌年以降も引き継ぎつつ、すべての変異発生数を推移行列で格納しておく。
ここで、ガンマ分布でさくっとやってしまうと、shape が発癌に必要なヒット回数、scale(rate) が1年あたりのヒット回数 になる。
# simulation n.ind <- 10000 # 人数 n.hit <- 10 # 発癌に必要なヒット回数 age.box <- 1:100 m2 <- 0.2 # 1年あたりのヒット回数。ポアソン分布の平均パラメータ first.cancer <- rep(-9,n.ind) for(i in 1:n.ind){ events <- rpois(length(age.box),m2) cumsum.events <- cumsum(events) tmp <- which(cumsum.events >= n.hit) if(length(tmp)>0){ first.cancer[i] <- tmp[1] - 1 } } # poisson model d <- dpois(0:(n.hit-1), m2) # less than n.hit probability d <- c(d, 1-sum(d)) # add more than n.hit probability # transition matrix # All transition patterns are considered. idx <- mapply(function(z) seq(z, n.hit+z, by=1), 0:n.hit) idx <- replace(idx, idx>n.hit, n.hit) dage <- 0:100 res <- matrix(0, 0:n.hit, nr=length(dage), nc=n.hit+1) res[1,] <- d for(i in 1:(length(dage)-1)){ trans <- outer(res[i,], d) #trans[lower.tri(trans)] <- 0 res[i+1,] <- tapply(c(trans), c(idx), sum) } par(mar=c(5, 5, 2.2, 0.5), cex.axis=1.5, cex.lab=1.5, font=2) hist(first.cancer,breaks=seq(from=-10,to=100,by=10), freq=FALSE, col=2, density=10, xlab="first cancer age") lines(head(dage, -1), diff(res[, n.hit+1]), col=4, lty=3, lwd=3) dg <- dgamma(dage, n.hit, rate=m2) lines(dage, dg, col=3, lwd=3) legend("topleft", legend=c("simulation", "pois model", "gamma model"), col=c(2,4,3), pch=15)
さて、ここで実際の癌の発生状況を確認してみる。国立がん研究センターから集計表をダウンロードして適当にネ申エクセルを使う部分だけtxt にする。
d <- read.delim("cancer.txt", check.names=FALSE, stringsAsFactors=FALSE) ye <- unique(d$診断年) idx <- grep("^C\\d+", d) m <- mapply(function(z) d[(idx[z]-2):(idx[z+1]-3)], head(seq(idx), -1)) ageidx <- grep("歳", colnames(d)) agename <- colnames(d)[ageidx] age <- sapply(sapply(str_extract_all(agename, "\\d+"), as.numeric), mean) dat <- mapply(function(z) split(z, z$性別), split(d, d$コード)) dat[[14]] <- list(dat[[14]][[1]], dat[[15]][[1]], dat[[16]][[1]]) dat[[13]] <- list(dat[[13]][[1]], dat[[17]][[1]], dat[[18]][[1]]) dat[15:16] <- dat[17:18] <- NULL yl <- range(0, d[,ageidx]/rowSums(d[,ageidx])) # 確率密度 cols <- jet.colors(nrow(tmp)) i5 <-ye %% 5 == 0 par(mfrow=c(length(dat), 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2, las=1) for(i in seq(dat)){ for(j in c(1,3,2)){ tmp <- dat[[i]][[j]] matplot(age, t(tmp[,ageidx]/rowSums(tmp[,ageidx])), type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3) title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2) legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2) } } # 発生数 par(mfrow=c(length(dat), 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2) for(i in seq(dat)){ #par(mfrow=c(1, 3)) yl <- range(0, do.call(rbind, mapply(function(z) z[,ageidx], dat[[i]], SIMPLIFY=FALSE))) for(j in c(1,3,2)){ tmp <- dat[[i]][[j]] matplot(age, t(tmp[,ageidx]), type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3) title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2) legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2) } }
さてここで、癌発生年齢の経験分布から、ガンマ分布のパラメータを適当に当てはめて推定してみる。ガンマ分布のパラメータは発癌までの回数と、1年あたりのヒット回数なので、これを比較すると癌になるまでに要する変異数的なものがわかる。
実際には、変異はDNA損傷修復機構で修復されるし、実際に発癌に寄与する変異(driver mutation) は固有の遺伝子のクリティカルな変異数個であることが多い(大腸がんや脳腫瘍など)。
若年での発症が極端には多くない大腸がんだけ取り出してみる。というのも、例えば白血病や脳・中枢神経系は10歳以下くらいの小児に後発するタイプの悪性腫瘍があるため、単純なガンマ分布だけでのフィッテイングだと当てはまりがよくない。また、若年発症するものは多くが遺伝が強く関係するタイプ(乳癌や子宮体癌のBRCA1/2 など)(ならば大腸がんだってFAP があるのでは、と思ったけど統計上は若年層はあまりいないようである)。皮膚がんや肝臓がんあたりは遺伝というより紫外線や肝炎ウイルスの寄与が大きいのでこれらはよりよくガンマ分布に当てはまりそうだと思ったけどそうでもなかった。
推定の初期値によるのかもしれないが、ヒット回数は12回から25回くらいが多いようである。ヒット速度は0.2 から0.3 くらい。
par(mfrow=c(1, 3), mar=c(2.5, 3, 2.2, 0.5), cex.axis=1.5, font=2, las=1) #for(i in seq(dat)){ for(j in c(1,3,2)){ tmp <- dat[[i]][[j]] y <- t(tmp[,ageidx]/rowSums(tmp[,ageidx])) s1 <- mapply(function(z){ ytmp <- y[,z] resid <- function(par){ yhat <- dgamma(age, par[1], par[2]) sum((ytmp-yhat)^2) } optim(init, resid, method="BFGS")$par}, seq(ncol(y))) dg <- mapply(function(z) dgamma(age, s1[1, z], s1[2, z]), seq(ye)) dg <- sweep(dg, 2, colSums(dg), "/") yl <- range(0, do.call(rbind, mapply(function(z) z[,ageidx]/rowSums(z[,ageidx]), dat[[i]], SIMPLIFY=FALSE))) tmp <- dat[[i]][[j]] hoge <- t(tmp[,ageidx]/rowSums(tmp[,ageidx])) matplot(age, dg, type="l", lty=1, col=cols, ylab="Density", ylim=yl, lwd=3) title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2) legend("topleft", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2) } } par(mar=c(5, 5, 2.2, 0.5), cex.axis=1.5, cex.lab=1.5, font=2) plot(t(s1), pch=15, col=cols, xlab="発癌までのヒット回数", ylab="1年あたりのヒット回数") legend("bottomright", legend=ye[i5], col=cols[i5], pch=15, ncol=2, cex=1.2) title(sprintf("%s: %s", tmp[1,2], tmp[1,4]), cex.main=2)
