ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2020年12月21日

対数ロジスティック分布をStanに実装する

 この記事は,Stan Advent Calendar 2020の21日目の記事です。この記事では,生存時間解析 (survival analysis; 継続時間解析とも) でよく使われる対数ロジスティック分布を使ったモデリングをStanで行う方法について解説します。ちなみに,一般化ガンマ分布の記事も以前書きました。この記事を読んで得られる学びとしては,(1)ロジスティック分布と対数ロジスティック分布の関係が分かる,(2)自作関数をStanで作る方法が分かる,(3)対数〇〇分布の代わりに対数変換したデータ+〇〇分布を使ってモデル化するとパラメータの推定はうまくいくがWAICの計算を間違えてしまうことがある,あたりが挙げられるかもしれません。

1. ロジスティック分布と対数ロジスティック分布

 対数ロジスティック分布は,ある確率変数Xの対数ln(X)がロジスティック分布に従うときに,元の確率変数Xが従う分布です。逆にいえば,ロジスティック分布に従う変数を指数変換した変数は対数ロジスティック分布に従います。正規分布と対数正規分布の関係と同じですね。ロジスティック分布の台 (実現値が取り得る範囲)は全ての実数ですが,対数ロジスティック分布の台は正の実数です。それゆえに,対数ロジスティック分布は非負の値となる生存時間や継続時間の確率分布の候補になる訳ですね。ロジスティック分布と対数ロジスティック分布の見た目はこんな感じです。

loglogidist.png

という訳で,まずはロジスティック分布の説明から。ロジスティック分布の確率密度関数は

f(x∣μ,σ)=exp⁡(-(x-μ)/σ)/(σ{1+exp⁡(-(x-μ)/σ) }^2 )

と表されます。μは位置パラメータ,σ (> 0) はスケールパラメータです。ちなみに,累積分布関数は

F(x∣μ,σ)=1/(1+exp⁡(-(x-μ)/σ) )

と書けます。この式に見覚えのある方も多いかと思いますが,μ = 0,σ = 1のとき (標準ロジスティック分布),累積分布関数は

F(x∣μ=0,σ=1)=1/(1+e^(-x) )=e^x/(1+e^x )

となりロジスティック関数と一致します。これがロジスティック分布という名前の由来ですね。標準ロジスティック分布とロジスティック関数の関係は,標準正規分布とその累積分布関数の関係と同じです。ロジスティック分布はStanではlogistic_lpdf()などの関数で実装されています (サンプリングステートメントを使って「Y ~ logistic(mu, sigma);」と書ける)。

 ある確率変数Xがロジスティック分布に従うとき,exp(X)は対数ロジスティック分布に従います。スケールパラメータα(> 0)と形状パラメータβ(> 0)を使って,対数ロジスティック分布の確率密度関数は

f(x∣α,β)=(β/α (x/α)^(β-1))/{1+(x/α)^β }^2

累積分布関数は

f(x∣α,β)=1/(1+(x/α)^(-β) )=x^β/(α^β+x^β )

と表されます。α = exp(μ),β = 1/σを代入すればμとσを使った表現もできますが,αとβを使った方がシンプルな式で表せます。どちらのパラメータ化を使うかは目的によって変えると良いでしょう。

2. 対数ロジスティック分布をStanに実装

 対数ロジスティック分布の関数はStanにはデフォルトで用意されていないので,(1) αとβでパラメータ化した対数ロジスティック分布を自分で定義する方法,(2) μとσでパラメータ化した対数ロジスティック分布を自分で定義する方法,(3) ロジスティック分布を利用する方法,の3通りの方法で実装してみます。3番目のロジスティック分布を使ったやり方は一番簡単ですが,後ほど示すように,WAICの計算などのために尤度を使う必要があるときには正しい結果を返さないので注意が必要です。結果を比較できるように,これらを1つのStanコードにまとめて書きます。WAICも計算します。Stanコードは以下の通り。

Stanコード: loglogistic.stan
functions{
  //1. αとβでパラメータ化
  real loglogistic1_lpdf(real x, real scale, real shape){
    real lpdf;
    lpdf = log(shape) - log(scale) + (shape-1)*log(x/scale) - 2*log(1+(x/scale)^shape);
    return(lpdf);
  }
  
  //2. μとσでパラメータ化
  real loglogistic2_lpdf(real x, real mu, real sigma){
    real lpdf;
    lpdf = -log(sigma) - mu + (1/sigma-1)*(log(x) - mu) - 2*log(1+(x/exp(mu))^(1/sigma));
    return(lpdf);
  }
}

data{
  int N;
  vector[N] Y1;
  vector[N] Y2;
  vector[N] Y3;
}

parameters{
  vector[3] mu;
  vector<lower=0>[3] sigma;
}

transformed parameters{
  vector<lower=0>[3] alpha;
  vector<lower=0>[3] beta;
  
  for(j in 1:3){
    alpha[j] = exp(mu[j]);
    beta[j] = 1/sigma[j];
  }
}

model{
  //prior
  mu ~ uniform(-10,10);
  sigma ~ uniform(0,10);
  
  //model
  for(i in 1:N){
    //1. αとβでパラメータ化
    Y1[i] ~ loglogistic1(alpha[1], beta[1]);
    
    //2. μとσでパラメータ化
    Y2[i] ~ loglogistic2(mu[2],sigma[2]);
    
    //3. ロジスティック分布を利用
    log(Y3[i]) ~ logistic(mu[3],sigma[3]);
  }
}

generated quantities{
  vector[N] log_lik_1;
  vector[N] log_lik_2;
  vector[N] log_lik_3;
  
  //WAICの計算のために対数尤度を計算する
  
  for(i in 1:N){
    //1. αとβでパラメータ化
    log_lik_1[i] = loglogistic1_lpdf(Y1[i] | alpha[1], beta[1]);
    
    //2. μとσでパラメータ化
    log_lik_2[i] = loglogistic2_lpdf(Y2[i] | mu[2],sigma[2]);
    
    //3. ロジスティック分布を利用 (間違った結果になります)
    log_lik_3[i] = logistic_lpdf(log(Y3[i]) | exp(mu[3]), 1/sigma[3]);
  }
}

 データとして与えたY1,Y2,Y3に対し,それぞれ(1)〜(3)の方法で対数ロジスティック分布を当てはめてパラメータを推定するコードです。4つの3次元ベクトルmu, sigma, alpha, betaの各要素にそれぞれの方法で推定されたパラメータの値が格納されます。functionsブロックでは,2通りのパラメータ化の方法で対数ロジスティック分布の対数確率密度関数 (lpdf) を定義しています (上述した確率密度関数の対数をとって整理したもの)。WAICを計算するために,generated quantitiesブロックで対数尤度の計算も行っています。ここではあえてmodelブロックと対応した書き方をしています (つまり,Y1〜Y3を使って3通りの方法で対数尤度を計算している)。Y1〜Y3を全て同じデータにしておけば,それぞれの方法の結果を比較することができます。

3. 架空のデータを使って走らせてみる

 それでは,Rで架空のデータを作って実際にパラメータを推定してみましょう。真値がμ=0.5 (α = 1.65),σ = 2.0 (β = 0.5) である対数ロジスティック分布からの乱数を5,000個発生させてみます。対数ロジスティック分布からの乱数を得るには,rlogis()で得たロジスティック分布からの乱数をexp()で指数変換するか,actuar::rllogis()などのパッケージに含まれる関数を使えばokです。Rコードと推定結果は以下の通り。

Rコード: loglogistic.R
library(tidyverse) #version 1.2.1
library(rstan) #version 2.19.2
library(loo) #2.3.1

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## 架空のデータを作る ----

N <- 5000
mu <- 0.5
sigma <- 2.0

set.seed(8823)

Y <- rlogis(N, location = mu, scale = sigma) %>% exp()

#もしくは
#library(actuar)
#Y <- rllogis(N, shape = 1/sigma, scale = exp(mu))

## Stanに渡す ----
sm <- stan_model("loglogistic.stan")

warmup <- 1000
iter <- 2500 + warmup
t0<-proc.time()
fit_loglogis <- sampling(sm,data=list(N=N,Y1=Y,Y2=Y,Y3=Y)
                ,seed=8823,chains=4,iter=iter,warmup=warmup
)

#saveRDS(fit_loglogis,"fit_loglogis.rds")

## 推定結果 ----
s_loglogis <- summary(fit_loglogis)$summary
round(s_loglogis,3)[1:12,]

#           mean se_mean    sd  2.5%   25%   50%   75% 97.5%    n_eff Rhat
# mu[1]    0.534   0.000 0.048 0.440 0.503 0.535 0.566 0.628 13459.11    1
# mu[2]    0.534   0.000 0.049 0.437 0.500 0.534 0.567 0.629 13529.81    1
# mu[3]    0.534   0.000 0.049 0.438 0.501 0.533 0.567 0.630 12189.30    1
# sigma[1] 1.991   0.000 0.023 1.946 1.975 1.990 2.006 2.036 13718.75    1
# sigma[2] 1.991   0.000 0.024 1.945 1.975 1.991 2.007 2.038 12643.81    1
# sigma[3] 1.991   0.000 0.023 1.945 1.975 1.990 2.006 2.037 14590.76    1
# alpha[1] 1.708   0.001 0.082 1.553 1.653 1.707 1.762 1.874 13440.38    1
# alpha[2] 1.707   0.001 0.084 1.548 1.649 1.706 1.763 1.877 13436.23    1
# alpha[3] 1.708   0.001 0.084 1.549 1.650 1.705 1.763 1.878 12220.63    1
# beta[1]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 13704.58    1
# beta[2]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 12601.31    1
# beta[3]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 14593.66    1

## WAICの計算
e_logliks <- rstan::extract(fit_loglogis, pars = c("log_lik_1", "log_lik_2", "log_lik_3"))
waic(e_logliks$log_lik_1) # WAIC = 32194.3 (SE = 527.7)
waic(e_logliks$log_lik_2) # WAIC = 32194.5 (SE = 527.8)
waic(e_logliks$log_lik_3) # WAIC = 53456.1 (SE = 662.6)

 パラメータの推定結果はどの方法でもほぼ同じ値でありかつ真値に近いので,うまく推定できていると見てよさそうです。一方,WAICの結果を見ると,方法(1)と方法(2)ではほとんど値が同じですが,ロジスティック分布を利用した方法(3)だけ結果が大きく異なっています。同じデータ・同じモデルなのにWAICが一致しないのはおかしいですね。これは,対数ロジスティック分布を使ったモデルではYの確率密度を使って尤度が計算されるのに対し,ロジスティック分布を使ったモデルではlog(Y)をデータとして与えたときの確率密度を使って尤度が計算されてしまうことに由来します。与えているデータが違うのだから,結果が変わるのは当たり前ですね (もっとシンプルな例として,データのスケールを変えると確率密度や尤度が変わる,ということを考えると分かりやすいかもしれません)。サンプリングするときは方法3のようにlog(Y[i]) ~ logistic(mu, sigma); という書き方をしても問題ありませんし,おそらくその方が自作関数を使うよりも計算が高速で収束も安定すると思われますが,尤度を正確に計算する必要があるときには書き方に気を付けましょう。これは対数正規分布を使うときも同様で,サンプリング時にはlog(Y) ~ normal(mu, sigma); と書いてもYが対数正規分布に従うことを表せますが,尤度を計算するときにはlog_lik[i] = lognormal(Y[i] | mu, sigma); という書き方をしなければなりません。Be careful and enjoy!

posted by mutopsy at 00:00 | 統計