【Stan】【農学】平均と標準偏差でガンマ分布を定義する：混合ガンマ分布による除草剤抵抗性雑草の割合推定

複数の確率分布を混合した分布（混合分布）で集団における各クラスタの割合推定をすることができます。

解析には２つのガンマ分布の混合分布を用いますが、その際、ガンマ分布のパラメータとして「平均・標準偏差」を用いるのが便利です。

今回は、農業雑草の研究をしている私のラボにおける仮想データから紹介します。

除草剤抵抗性雑草について

農業の現場において、作物の収量低下をもたらす雑草は、現在では除草剤処理により抑制することが多いです。

しかし、一つの除草剤の繰り返しの使用や、規定の薬量濃度に達しない低薬量散布などの要因により、 「除草剤抵抗性」形質を持つ個体が出現し、世界中で大きな農業被害をもたらしています。

ある個体が除草剤抵抗性の形質を有しているかを室内実験から判断するには、

・除草剤抵抗性を有しないことが既知の（感受性の）対照群
・ある水田圃場から採取した対象の個体群

のそれぞれに除草剤処理を施し、 t検定等を行い判断します。

すなわち、「抵抗性 / 感受性」の評価を一個体単位で行うのではなく、個体数に反復を取った水田圃場等の集団単位で行います。

これは、たとえばある個体が抵抗性（＝抵抗性遺伝子をいくらか有している個体）であったとしても、
薬量処理濃度や、測定する評価指標（葉・根・重量etc.）によっては感受性と判断がつかないほどに生育不良を起こすことも多いためです。

データの分布

本題です。

今回は、上記のような個体群単位の情報が得られなかった場合を想定します。

たとえば、海外からの輸入穀物に雑草の種子がよく混入しているのですが、
これに含まれる除草剤抵抗性個体の割合を推定したい、とします。

ここにもはや圃場単位の情報は存在しないため、一個体ごとに抵抗性か感受性かを判断しなければなりません。

これらの輸入穀物混入種子を発芽させ、除草剤処理を施しました。
生育の評価指標として、「薬剤処理後の根の伸長量」を用いました。

200個体に対して実験を施し、得られた200行のデータの分布は以下のようになりました。

根の伸長量が小さい除草剤感受性型と、根の伸長量が大きい除草剤抵抗性型が混在することによる双峰型の分布になりました。

これを２つのガンマ分布の混合分布と見立て、除草剤抵抗性集団の割合を推定します。

なお、シミュレーションデータを生成した２つのガンマ分布の真値は以下のようにしました（コードはページ最後に載せています）。

パラメータはshapeパラメータとrateパラメータで指定しており、 図の二つの分布を、除草剤抵抗性個体（図の右の山）の割合aを0.7として混合しています。

モデルの記述

data{
  int N;          // サンプルサイズ
  vector[N] Y;    // 根の伸長量
}

parameters{
  real<lower=0, upper=1> a;   // 抵抗性個体の割合
  positive_ordered[2] shape;  // shapeパラメータ
  vector<lower=0>[2] rate;    // rateパラメータ
}

model{
  // 対数尤度の計算
  for(n in 1:N)
    target += log_sum_exp(
      log1m(a)   + gamma_lpdf(Y[n] | shape[1], rate[1]),  // 感受性個体
      log(a) + gamma_lpdf(Y[n] | shape[2], rate[2])       // 抵抗性個体
    );
}

感受性（図の左の山）のガンマ分布のパラメータを(shape1, rate1)、
抵抗性（図の右の山）を(shape2, rate2)とし、
抵抗性の割合をaとしたとき、

ある一サンプルのデータが得られる尤度pが

であり、これの対数尤度が

ですので、
log_sum_exp関数の定義である

　　log_sum_exp(x, y) = log (exp(x) + exp(y))

に従いmodelブロックで記述し、対数尤度を足し上げています。
なお、log1m(a) = log(1-a) です。

ラベルスイッチングを避けるため、正のパラメータであるshapeパラメータに shape[1] < shape[2] となる制約をpositive_ordered型で宣言することにより課しています。

これをRからキックすればよいわけですが、結果から言うと収束しません。

トレースプロットを見ると、以下のようになっています。

rate1が彼方に飛んでいってしまっているようです。

seed値や元データによってはshapeパラメータが彼方に飛んでいったりもします。

こういったあり得ない数値に飛ぶことを防止するため、
shapeパラメータとrateパラメータに範囲の制約を課せばいいのですが、
ガンマ分布のパラメータの解釈は少しわかりづらいという問題があります。

スクリプトを他者に見せる際にも、この値はどこから出てきたのかな？と思われそうでもあります。

そこで、より直感的なパラメータ、
すなわち平均と標準偏差を用いてガンマ分布を定義し直します。

平均と標準偏差でガンマ分布を定義する

shapeパラメータ、rateパラメータ のガンマ分布について、
平均muと標準偏差sigmaはそれぞれ
　 　
で表されます。

これを変形すると、shapeパラメータとrateパラメータは、平均と標準偏差を用いて
　 　
となります。

これを用いてモデルを記述し、各パラメータに制約を入れます。

data{
  int N;          // サンプルサイズ
  vector[N] Y;    // 根の伸長量
}

parameters{
  real<lower=0, upper=1> a;           // 抵抗性個体の割合
  real<lower=0, upper=15> mu1;        // 感受性個体の平均
  real<lower=0, upper=100> mu2;       // 抵抗性個体の平均
  real<lower=0, upper=15> sigma1;     // 感受性個体の標準偏差
  real<lower=0, upper=100> sigma2;    // 抵抗性個体の標準偏差
}

transformed parameters{
  real shape1 = square(mu1 / sigma1);  // 感受性個体のshapeパラメータ
  real shape2 = square(mu2 / sigma2);  // 抵抗性個体のshapeパラメータ
  real rate1 = mu1 / square(sigma1);   // 感受性個体のrateパラメータ
  real rate2 = mu2 / square(sigma2);   // 抵抗性個体のrateパラメータ
}

model{
  // 対数尤度の計算
  for(n in 1:N)
    target += log_sum_exp(
      log1m(a)   + gamma_lpdf(Y[n] | shape1, rate1),
      log(a) + gamma_lpdf(Y[n] | shape2, rate2)
    );
}  

パラメータの制約なしではやはり収束しませんが、
parametersブロックに示したような制約を入れることで収束します。

ここでは、
「感受性個体の根の伸長量の平均は15mm以下であり、抵抗性個体は100mm以下、 またそれぞれのバラツキもそれより大きくなることはない」
という意味の制約を課しています。

shapeパラメータとrateパラメータを直接記述するよりも、
平均と標準偏差を用いた定義の仕方では、
より直感的な解釈が可能になっていると分かります。

今回は単にモデルを収束させるための範囲制限でしたが、
たとえば平均と標準偏差を用いて事前分布を記述することも容易になります。

ちなみに、推定結果の事後分布は以下のようになります。

すべてのパラメータはその95%区間に真値を含んでおり、良い精度で推定できています。
除草剤抵抗性個体の推定割合aは、事後分布の中央値で69%であると評価されました。

参考

・StanとRでベイズ統計モデリング　松浦健太郎
　混合分布やlog_sum_exp関数について詳しく説明されています。

Rコード

最後に、使用したRコードを載せておきます。

library(rstan)
library(ggplot2)
set.seed(100)

# パラメータ真値の設定
a <- 0.7
N <- 200
shape1 <- 2
shape2 <- 5
rate1 <- 0.7
rate2 <- 0.2

 # 混合前の分布図示
curve(dgamma(x, shape1, rate1), 0, 60, ylab="", xlab="")
curve(dgamma(x, shape2, rate2), 0, 60, add=T)

y1 <- rgamma(N*(1-a), shape1, rate1)
y2 <- rgamma(N*a, shape2, rate2)

d <- data.frame(Y=c(y1, y2)) # シミュレーションデータセット

# データの把握
ggplot(d, aes(x=Y, y=..density..))+
  theme_test()+
  geom_histogram(fill="red", alpha=0.5)

# Stanのキック
data <- list(N=nrow(d), Y=d$Y)

# shapeパラメータ、rateパラメータによる記述
fit1 <- stan(file="params_default.stan", data=data, seed=123)
stan_trace(fit1)  # 収束しない

# 平均と標準偏差からパラメータを指定
fit2 <- stan(file="params_mean_sd.stan", data=data, seed=123)
stan_trace(fit2)
stan_hist(fit2)