多重共線性とか潜在変数とか共分散構造分析とか(2)

承前。モデルを見直してみた。

sem.bug

var
  N,    # number of data
  x1[N], x2[N], x3[N], x4[N], x5[N], x6[N], # data
  xi1[N], eta1[N],  # latent variables
  lambda11, lambda21, lambda31, lambda41,   # parameters
  kappa51, kappa61, mu,
  mu1[N], mu2[N], mu3[N], mu4[N], mu5[N], mu6[N],
  tau1, tau2, tau3, tau4, tau5, tau6,
  sigma1, sigma2, sigma3, sigma4, sigma5, sigma6

model {
  for (i in 1:N) {
    x1[i] ~ dnorm(mu1[i], tau1);
    mu1[i] <- lambda11 * xi1[i];
    x2[i] ~ dnorm(mu2[i], tau2);
    mu2[i] <- lambda21 * xi1[i];
    x3[i] ~ dnorm(mu3[i], tau3);
    mu3[i] <- lambda31 * xi1[i];
    x4[i] ~ dnorm(mu4[i], tau4);
    mu4[i] <- lambda41 * xi1[i];
    x5[i] ~ dnorm(mu5[i], tau5);
    mu5[i] <- kappa51 * eta1[i];
    x6[i] ~ dnorm(mu6[i], tau6);
    mu6[i] <- kappa61 * eta1[i];
    eta1[i] ~ dnorm(xi1[i], 1);
    xi1[i] ~ dnorm(mu, 1);
  }
  lambda11 ~ dgamma(1.0E-3, 1.0E-3);
  lambda21 ~ dgamma(1.0E-3, 1.0E-3);
  lambda31 ~ dgamma(1.0E-3, 1.0E-3);
  lambda41 ~ dgamma(1.0E-3, 1.0E-3);
  kappa51 ~ dgamma(1.0E-3, 1.0E-3);
  kappa61 ~ dgamma(1.0E-3, 1.0E-3);
  mu ~ dnorm(0, 1.0E-6);
  tau1 ~ dgamma(1.0E-3, 1.0E-3);
  tau2 ~ dgamma(1.0E-3, 1.0E-3);
  tau3 ~ dgamma(1.0E-3, 1.0E-3);
  tau4 ~ dgamma(1.0E-3, 1.0E-3);
  tau5 ~ dgamma(1.0E-3, 1.0E-3);
  tau6 ~ dgamma(1.0E-3, 1.0E-3);
  sigma1 <- 1/sqrt(tau1);
  sigma2 <- 1/sqrt(tau2);
  sigma3 <- 1/sqrt(tau3);
  sigma4 <- 1/sqrt(tau4);
  sigma5 <- 1/sqrt(tau5);
  sigma6 <- 1/sqrt(tau6);
}

パラメータの事前分布をガンマ分布とすることで正値に固定。一方、muは正規分布に。これで初期値の指定が不要になった。

library(rjags)
mcmc.data <- list(N = N, x1 = x1, x2 = x2, x3 = x3,
             x4 = x4, x5 = x5, x6 = x6)
mcmc.model <- jags.model("sem.bug", mcmc.data,
                         nchain = 3,
                         n.adapt = 10000)
vars <- c("lambda11", "lambda21", "lambda31",
          "lambda41", "kappa51", "kappa61",
          "sigma1", "sigma2", "sigma3", "sigma4",
          "sigma5", "sigma6", "mu")
mcmc.samples <- coda.samples(mcmc.model, vars,
                             n.iter = 20000, thin = 20)

結果

Iterations = 10020:30000
Thinning interval = 20 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean      SD  Naive SE Time-series SE
kappa51  2.01337 0.20272 0.0037011      0.0125964
kappa61  3.96202 0.38696 0.0070648      0.0254412
lambda11 0.70533 0.06788 0.0012392      0.0035674
lambda21 0.71944 0.07728 0.0014109      0.0036421
lambda31 1.42915 0.14336 0.0026174      0.0073583
lambda41 1.45978 0.15648 0.0028569      0.0073220
mu       0.05931 0.16088 0.0029373      0.0033445
sigma1   0.10931 0.04148 0.0007573      0.0011649
sigma2   0.29772 0.03956 0.0007223      0.0007045
sigma3   0.35253 0.07052 0.0012876      0.0017547
sigma4   0.56326 0.08367 0.0015276      0.0019117
sigma5   0.69995 0.34757 0.0063458      0.0244314
sigma6   0.87967 0.72036 0.0131520      0.0522985

2. Quantiles for each variable:

             2.5%      25%     50%    75%  97.5%
kappa51   1.63720  1.87076 1.99980 2.1521 2.4291
kappa61   3.25292  3.68748 3.94266 4.2239 4.7374
lambda11  0.58796  0.65860 0.69953 0.7465 0.8601
lambda21  0.58781  0.66598 0.71244 0.7653 0.8833
lambda31  1.18391  1.32766 1.41589 1.5188 1.7523
lambda41  1.18920  1.35072 1.44958 1.5521 1.7988
mu       -0.25829 -0.05286 0.06216 0.1692 0.3695
sigma1    0.03231  0.07888 0.11044 0.1377 0.1887
sigma2    0.23146  0.26966 0.29385 0.3216 0.3829
sigma3    0.20817  0.30843 0.35293 0.3977 0.4875
sigma4    0.41712  0.50489 0.55680 0.6141 0.7446
sigma5    0.04025  0.41204 0.83364 0.9566 1.1568
sigma6    0.03894  0.18897 0.67682 1.5982 2.1305

タグ：SEM MCMC