SEMの多重指標モデルとMIMICモデル

豊田先生らの「原因をさぐる統計学」の読了記念に記事にします。データ及びRの{lavaan}を使ったスクリプトは豊田先生の研究室のホームページからダウンロードできます。1992年の本なのにわざわざ2013年に{lavaan}でやった例を追加していて敬服します。今回は書籍のp.107の「食物摂取とがん」データの分析（柳井・高木 (1986)）についてやってみました。本持っていない人には分かりにくくてすみません。

以下、今回のネタ本と実装にあたって参考にした本です。

原因をさぐる統計学―共分散構造分析入門 (ブルーバックス)

MIMICモデル

先にStanで組みやすいMIMICモデルから。MIMICモデルは矢印の流れが左から右へと一方向になっており、予測に適したモデルとなっているとのこと。今回の例では、摂取カロリーや摂取肉類などの説明変数（Y[1]～Y[4]）から「下部消化管のガン傾向」という潜在変数が回帰で決まるとし、そこからまた回帰で大腸ガン（Y[5]）、直腸ガン（Y[6]）の各死亡率が回帰で決まると考えます。Stanでの実装は以下になりました。

data {
   int<lower=1> N;
   vector[N] Y[6];
}

parameters {
   vector[N] mu;
   real beta[5];
   real<lower=0> sigma[3];
}

model {
   mu ~ normal(beta[1]*Y[1] + beta[2]*Y[2] + beta[3]*Y[3] + beta[4]*Y[4], sigma[1]);
   Y[5] ~ normal(      1*mu, sigma[2]);
   Y[6] ~ normal(beta[5]*mu, sigma[3]);
}

• 14行目: 出力の片方の係数は1にします。ここもbeta[*]をかけてしまうと、「muが定数倍する」のと「beta[*]とbeta[5]が定数倍する」のが等価になってしまい、決まらなくなってしまうためです。これはmuのスケールをY[5]にあわせてる意味合いです。なお、vector型を積極的に使っているのでfor文が省略されていても要素ごとに評価が行われています。少しコードが簡潔になる分、分かりにくくなっています。ご留意下さい。

キックするRコードは以下です。

library(rstan)
stan_fit <- stan(file='MIMIC.stan', chains=0)
cancer.dat <- read.csv('input/Cancer_raw.csv')
N <- nrow(cancer.dat)
Y <- t(scale(cancer.dat))
data <- list(N=N, Y=Y)
fit <- stan(fit=stan_fit, data=data, iter=2500, warmup=500, thin=1, chains=3, seed=123)

• 5行目では入力となるデータを標準化しています。特に中央化することで切片項を削ることは常套手段のようです。 計算時間は1chainあたり5秒以内でした。回帰係数の信用区間が求めることができたりしてハッピーです。

多重指標モデル

多重指標モデルは「洋食傾向」という潜在変数を考えて、そこから摂取カロリーや摂取肉類など（Y[1]～Y[4]）が1つ1つ回帰で決まると考えます。また「洋食傾向」から「下部消化管のガン傾向」という潜在変数が回帰で決まるとし、あとはMIMICと同様です。重回帰がなく単回帰をつなげたものになっており、解釈が容易となるとのこと。実装は以下になりました。

data {
  int<lower=1> N;
  vector[6] Y[N];
}

parameters {
  vector[N] mu1;
  vector[N] mu2;
  vector[6] a;
  real<lower=0> gam;
  vector<lower=0>[6] lam;
  vector<lower=0>[6] sigma;
}

model {
  mu1 ~ normal(0, 1);
  mu2 ~ normal(gam*mu1, 1);
  lam ~ normal(1, 1);
  for (n in 1:N) {
    Y[n,1] ~ normal(a[1] + lam[1]*mu1[n], sigma[1]);
    Y[n,2] ~ normal(a[2] + lam[2]*mu1[n], sigma[2]);
    Y[n,3] ~ normal(a[3] + lam[3]*mu1[n], sigma[3]);
    Y[n,4] ~ normal(a[4] + lam[4]*mu1[n], sigma[4]);
    Y[n,5] ~ normal(a[5] + lam[5]*mu2[n], sigma[5]);
    Y[n,6] ~ normal(a[6] + lam[6]*mu2[n], sigma[6]);
  }
}

lamを正の値に制限しているのは, mu1とlamが同じ絶対値のまま両方の符号を反転させた値も解になるためです。mu1とlamの事前分布は以下の本を参照しました。

Applied Bayesian Modelling (Wiley Series in Probability and Statistics Book 595) (English Edition)

分散共分散行列を使ったモデル

メモ。またいつか更新します。

data {
  int<lower=1> N;
  vector[6] Y[N];
}

transformed data {
  vector[2] Zero;
  Zero = rep_vector(0, 2);
}

parameters {
  vector[2] mu[N];
  real<lower=-1, upper=1> rho;
  vector[6] a;
  vector<lower=0>[6] lam;
  vector<lower=0>[6] sigma;
}

transformed parameters {
  matrix[2,2] cov;
  cov[1,1] = 1; cov[1,2] = rho;
  cov[2,1] = rho; cov[2,2] = 1;
}

model {
  mu ~ multi_normal(Zero, cov);
  lam ~ normal(1, 1);
  for (n in 1:N) {
    Y[n,1] ~ normal(a[1] + lam[1]*mu[n,1], sigma[1]);
    Y[n,2] ~ normal(a[2] + lam[2]*mu[n,1], sigma[2]);
    Y[n,3] ~ normal(a[3] + lam[3]*mu[n,1], sigma[3]);
    Y[n,4] ~ normal(a[4] + lam[4]*mu[n,1], sigma[4]);
    Y[n,5] ~ normal(a[5] + lam[5]*mu[n,2], sigma[5]);
    Y[n,6] ~ normal(a[6] + lam[6]*mu[n,2], sigma[6]);
  }
}

多重指標モデルと分散共分散行列を使ったモデルをkickするRコードは以下の通り。

library(rstan)

# stanmodel <- stan_model(file='model/model_MI.stan')
stanmodel <- stan_model(file='model/model_cov.stan')

cancer.dat <- read.csv('input/Cancer_raw.csv')
N <- nrow(cancer.dat)
Y <- scale(cancer.dat, center=FALSE)
data <- list(N=N, Y=Y)
fit <- sampling(stanmodel, data=data, seed=123)