順序制約のあるrstan

推定するパラメータに順序制約があるとき、rstan ではordered、正に限定するならpositive_ordered が使える。
　
例えばYaegashi (J Hum Genet. 1998;43(2):85-90.)らは、出生時の染色体異常が母体年齢に応じてどう推移するかをデータを取って調べている。
ここで、筆者らは というモデルを考えている。

txt <- "
   cases +21
35   736   1
36   850   4
37   824   8
38   800   2
39   713   8
40   626   6
41   417   3
42   255   6
43   146   2
44    87   2
45    22   0
46     7   0
47     1   0
"
yaegashi <- read.table(text=txt, check.names=FALSE)
y <- head(yaegashi$"+21"/yaegashi$cases, -3)
x <- 35:44
lm0 <- lm(log(y) ~ x)

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, cex.axis=1.5, las=1)
plot(x, y, log="y", ylim=c(0.5, 100)/1000, yaxt="n", pch=16, xlab="Maternal age", ylab="Rate per 1000", cex=2)
axis(2, at=c(1,5,10,50,100)/1000, labels=c(1,5,10,50,100))
lines(x, exp(cbind(1, x) %*% lm0$coefficients), lwd=3)
# exp(lm0$fitted.value)

lm0

Call:
lm(formula = log(y) ~ x)

Coefficients:
(Intercept)            x  
   -14.5127       0.2447  

なので、係数 は論文のTrisomy 21 の値 0.245 と一致している。
図も回帰線とプロットをそれなりに真似るとこんな感じ。

　
ここで、45週以降のTrisomy 21 の発生は0件なので、 の回帰では45週以降の発生率の観測値が0 になるので、線形回帰ではよろしくない。というわけで45週以降は欠測値として回帰したが、rstan による状態空間モデルにより、45週以降の発生率も推定してみる。
　
芸がなく2階差分でモデルを作る。
母体年齢 に応じた発生率は で表されるとして、ある年齢の妊婦 人から実際に観測される症例数は二項分布に従う、とすれば、以下のように書ける。
（ここではまだ順序制約はいれていない）

code <- "
data{
  int Nage;
  int age[Nage];
  int Pos[Nage];
  int Ntest[Nage];
}
parameters{
  real<lower=0, upper=1> p[Nage];
  real<lower=0> s;
}
model{
  for(i in 3:Nage){
    p[i] ~ normal(2*p[i-1] + p[i-2], s);
  }
  for(i in 1:Nage){
    Pos[i] ~ binomial(Ntest[i], p[i]);
  }
}

"

# sampling
standata <- list(Nage=nrow(yaegashi), age=as.numeric(rownames(yaegashi)), Pos=yaegashi$"+21", Ntest=yaegashi$cases)
fit_yaegashi <- sampling(m, standata, warmup=1000, iter=5000, seed=1234)

　
推定結果は下に図示するが、年齢に応じて凸凹している。というのも、38歳で突如、観測された発生率が落ち込むので、これに引きずられて推定値も下がる。
というわけで、年齢に応じて発生率が凸凹するのはなんか気持ち悪い、と思ったときに、年齢に応じた順序制約を入れる。
これはrstan ではordered もしくはpositive_ordered で作ることができる。発生率 は0以上1以下の確率なので、positive_ordered が良さそう、と思って

positive_ordered<lower=0, upper=1>

みたいなことをするのだが、rstan はこの書き方を受け付けないらしい。
とういわけでordered からtransformed parameters で変換するのが手らしい。しかしコンパイル時に警告をはく。無視。
　
実数値を[0,1] に変換したければ、ロジットの出番である。

# positive ordered に書き換えたコード
code <- "
data{
  int Nage;
  int age[Nage];
  int Pos[Nage];
  int Ntest[Nage];
}
parameters{
  //real<lower=0, upper=1> p[Nage];
  ordered[Nage] pp;
  real<lower=0> s;
}
transformed parameters{
  real<lower=0, upper=1> p[Nage];  // 変換されたパラメータ
  for(i in 1:Nage){
    p[i] = 1/(1+exp(-pp[i]));  // ロジットにすれば[0,1] に強制的に収まる
  }
}
model{
    for(i in 3:Nage){
      p[i] ~ normal(2*p[i-1] + p[i-2], s);
  }
  for(i in 1:Nage){
      Pos[i] ~ binomial(Ntest[i], p[i]);
    }
}

　
順序制約のないときはジグザグ凸凹な推定だったが、順序制約をいれると年齢に応じて発生率が上昇するわりときれいめな推定結果になった。
また、実際の観測数が0 である45週以上の推定も何とかできる。ただし、ここはということを利用しているだけなので、数が少ないため信用区間はかなり広い。