成長速度の変化点 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-05-31

読んだ Behav Ecol Sociobiol. 2007 Dec;62(2):245-253
dimorphic allometry という二形性の成長を見ていて、その変化点らしきところを推定する。
論文ではT. dichotomus septentrionalis と呼ばれる、いわゆるカブトムシの体の大きさ(body) と角(horn) についてデータをとっている。
　
データについて、対数をとったbody の大きさ $X$ と対数をとったhorn の長さ $Y$ について
$Y=\alpha_0+\alpha_1X+\alpha_2X^2+\epsilon$
という線形回帰(パラメータに関して言えば、線形)を考え、 $\alpha_2$ が有意に0から離れていれば、 $x^2$ の項がきいてくるので以降の解析は非線形でやります、といっているが、言いたいことは確かにわかるんだが、これってただ $x^n\hspace{3}(n\geq 2)$ の高次の項を含みうりますって言ってるだけで、パラメータ的に非線形であるかは否定できてないような気がしないでもない。
　
本筋ではないので解析にいくと、Fig 1 を見ると、なんとなくふたつの集団に分かれそうな点である。

このデータを真似たシミュレーションだと、kmeans で終わってしまったのだが、この集団を分けそうな $Y$ の値を回帰、というかstanで求める。
いま、 $X^0$ をswitch point として、ここを境にbody が影響してくるとすると
$Y=\beta_0+\beta_1X+\beta__2(X-X^0)D+\beta__3D+\epsilon$
ここで、 $D=\{\begin{matrix}1&if\hspace{3}X-X^0>0\\0&otherwise\end{matrix}$ である。
$\beta_2$ が有意に0から離れていれば、 $X^0$ を境に $X$ と $Y$ の関係式の傾きが変化する。
$\beta_3$ が有意に0から離れていれば、 $X^0$ を境に $Y$ は不連続になる。
$\beta_2$ か $\beta_3$ か両方が有意であれば、統計学的にdimorphic と言いたい。
　
$X^0$ についてやったので、 $Y^0$ も同様に
$X=\beta_0+\beta_1Y+\beta__2(Y-Y^0)D+\beta__3D+\epsilon$
として推定する。
$Y^0$ については、これより大きいhorn を持つ個体をMajor, 小さいhorn を持つ個体をMinor と定義する。
最適な $X^0$ と $Y^0$ は、 $X^0$ と $Y^0$ を動かしまくって決定係数 $R^2$ が最大となるものを採用するらしい。めっちゃ多重検定感()
　
Major とMinor な個体が決定したら、その群がallometry (相対的にめっちゃ成長する)かisometry (成長比率は普通)なのかを推定する。allometry 係数 $\alpha$ は、適当な係数 $a$ を用いて
$y=ax^{\alpha}$
$log{y}=log{a}+\alpha log{x}$
を回帰して推定する。
$\alpha>1$ ならば、 $y$ (horn) は $x$ (body) に比してクッソ大きい(positive allometry)
$a=1$ ならば、 $y$ と $x$ の関係は普通(isometry)
$0<\alpha<1$ ならば、 $y$ (horn) は $x$ (body) に比してクッソ小さい(negative allometry)
となる。
　
データを作る。

f <- function(x, a, b) a*x + b
a <- c(2.13, 0.79)
b <- c(-1.502, 0.369)

xl <- list(c(1.15, 1.4), c(1.3, 1.45))

n <- c(110, 491)
sd <- c(0.05, 0.03)
i <- 1
l <- mapply(function(i){
  x <- runif(n[i], min=xl[[i]][1], max=xl[[i]][2])
  hoge <- mapply(f, a=a[i], b=b[i], x=x) + rnorm(n[i], sd=sd[i])
  list(body=x, horn=hoge)
}, seq(n), SIMPLIFY=FALSE)

xlim <- c(0.9, 1.6)
ylim <- c(0.9, 1.6)
delta <- 0.03
cols <- c("lightblue", "green3")
plot(0, type="n", xlim=xlim, ylim=ylim, log="", xlab="Log body size", ylab="Log horn size")
for(i in seq(l)){
  x <- seq(xl[[i]][1]-delta, xl[[i]][2]+delta, length=100)
  y <- mapply(f, a=a[i], b=b[i], x=x)
  points(l[[i]]$body, l[[i]]$horn, col=cols[i], pch=16)
  lines(x, y)
}
abline(0, 1, lty=3)
legend("bottomright", legend=paste(c("Minor", "Major"), "Males"), pch=16, col=cols)

xy <- list(x=c(1.25, 1), y=c(1, 1.5))
for(i in 1:2) text(xy$x[i], xy$y[i], paste(c("Minor ","Major ")[i], "males:\ny=", a[i], "x+", b[i], sep=""), pos=4)

図はこんな感じ。

　
答えがわかっていない図はこんな感じ。これを解析しよう。

　
stan を使ってみる。 $D$ はstep 関数を使って書ける。if_else でゴリ押しもできる。

# morphometric.stan
data {
  int<lower=0>  n; # number of samples
  real<lower=0> Y[n]; # Y variable
  real<lower=0> X[n]; # X variable
}
parameters{
  real<lower=min(X), upper=max(X)> X0;   # switch point
  real b[4]; # cofficients b
  real<lower=0> s;
}
model {
  for(i in 1:n)
    Y[i] ~ normal(b[1] + b[2]*X[i] + (b[3]*(X[i] - X0) + b[4])*step(X[i] - X0), s);
}

　
それらしい $X^0$ が推定できた。しかし、 $\hat{R}$ は2 程度と収束が悪い。traceplot を見ると、chain 3 は良さそうでそれっぽいが、ほかはバクハツしている。

dat <- do.call(rbind.data.frame, l10)
xlim <- c(0, 30)
ylim <- c(0, 40)
plot(dat, xlim=xlim, ylim=ylim, pch=16, xlab="body size (mm)", ylab="horn size (mm)")

library(rstan)
stanmodel01 <- stan_model("morphometric.stan")
standata <- list(n = nrow(dat), X = dat$body, Y = dat$horn)

fit0 <- sampling(stanmodel01, data=standata, chains=3, warmup=500, iter=1500, seed=1234)
ex0 <- extract(fit0)

X0 <- median(ex0$X0)
b <- apply(ex0$b, 2, median)
b_quantile <- apply(ex0$b, 2, quantile, c(0.025, 0.975))

# 推定値のプロット用
x <- seq(min(dat$horn), max(dat$horn), length=1000)
y0 <- colSums(rbind(1, x, (x-X0)*0, 0)*b)
y1 <- colSums(rbind(1, x, (x-X0)*1, 1)*b)

w <- dat
plot(w, pch=16, xlim=xlim, ylim=ylim)
lines(x, y0)
lines(x, y1)
abline(v=X0, lty=3)
text(X0, par()$usr[3], substitute(italic(X^0)), pos=1, xpd=TRUE)
legend("bottomright", legend=substitute(italic(X^0)==x, list(x=round(X0, 2))))

　
$Y^0$ についても推定する。sampling 関数に渡すstandata を入れ替えるだけ。
これもいい感じの $Y^0$ は推定できたが、 $\hat{R}$ は6 とさらにひどい。

standata <- list(n = nrow(dat), Y = dat$body, X = dat$horn) # X Y 入れ替え
fit1 <- sampling(stanmodel01, data=standata, chains=3, warmup=500, iter=1500, seed=123)

ex1 <- extract(fit1)

Y0 <- median(ex1$X0)
b <- apply(ex0$b, 2, median)

plot(w, pch=15)
lines(x, y0)
lines(x, y1)
abline(v=X0, h=Y0, lty=3)
text(X0, par()$usr[3], substitute(italic(X^0)), pos=1, xpd=TRUE)
text(par()$usr[1], Y0, substitute(italic(Y^0)), pos=2, xpd=TRUE)
legend("bottomright", legend=substitute(atop(italic(X^0)==x, italic(Y^0)==y), list(x=round(X0, 2), y=round(Y0, 2))), adj=c(0, 0.2))

# Y0 classification
idx <- dat$horn > Y0
idx0 <- ifelse(idx, "major", "minor")

xlim <- c(0.9, 1.6)
ylim <- c(0.9, 1.6)
plot(w, log="xy", pch=idx+15, col=cols[idx+1], xlim=10^xlim, ylim=10^ylim)
abline(l[[1]]$coefficients)
abline(l[[2]]$coefficients)
abline(v=X0, h=Y0, lty=3)
text(X0, 10^par()$usr[3], substitute(italic(X^0)), pos=1, xpd=TRUE)
text(10^par()$usr[1], Y0, substitute(italic(Y^0)), pos=2, xpd=TRUE)
legend("bottomright", legend=names(pv), pch=15+1:0, col=cols[idx+1])

　
$\alpha$ について、Minor は2, Major は0.8 程度でどちらも有意になるので、Minor はpositive allometry, Major はnegative allometry となる。

alpha <- mapply(function(v){
  lm(formula = horn ~ body, data = log(w[v, ], 10))
}, split(seq(idx), idx0), SIMPLIFY=FALSE)
pv <- mapply(function(v) summary(v)$coefficients[2, "Pr(>|t|)"], alpha)

$major

Call:
lm(formula = horn ~ body, data = log(w[v, ], 10))

Coefficients:
(Intercept)         body  
     0.3086       0.8328  


$minor

Call:
lm(formula = horn ~ body, data = log(w[v, ], 10))

Coefficients:
(Intercept)         body  
     -1.261        1.934