Solved – How to specify the Wishart distribution scale matrix

jagspriorrwishart-distribution

I'm running the below Bayesian mixing model in R using the rjags package, but I am having difficultly in specifying the scale matrix for the Wishart distribution. Essentially, I want Sigma.inv to be a covariance matrix that accounts for instrumental measurement error as determined from repeat analysis of a control sample. I have chosen the Wishart distribution as the Sigma.inv prior, with the scale matrix taking the form of the covariance matrix of the control samples.However, this results in the model failing to converge, even after 1 million runs with 500,000 burn-in and 125 thinning length, and tiny errors are generated for the estimates of p.

So my question is should the scale matrix be in a different form? Or is there perhaps an alternative to the Wishart distribution that would be more appropriate in this case? I should add that the model works quite well if a use the identity matrix to scale the Wishart distribution, but of course this doesn't represent the measurement error that I want to incorporate.

Any help will be most appreciated.

model
{
 # Likelihood
for(i in 1:N) {
Y[i,1:J] ~ dmnorm(mu[i,],Sigma.inv)
mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]
}

# Priors on s
for(i in 1:N) {
for(k in 1:K) {
  s[k,1:J,i] ~ dmnorm(mu.s[,k],Sigma.s.inv[,,k])
  }
}

# Prior on p
for(i in 1:N) {
p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])
phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV
for(k in 1:K) {
  expphiV[i,k] <- exp(phiV[i,k])
}
}

# Prior on phi
for(i in 1:N) {
  for(k in 1:(K-1)) {
  phi[i,k] ~ dnorm(0,kappa.inv[k])
}
philast[i] <- -sum(phi[i,1:(K-1)])
}

for(k in 1:(K-1)) {
kappa.inv[k] ~ dgamma(2,1)
kappa[k] <- 1/kappa.inv[k]
}

# Prior on Sigma
Sigma.inv ~ dwish(R.Sigma,k.Sigma)
}

#

R Code to run model

#

# Y values
Y<-matrix(c(5.8565, 18.5979,    0.0045, 5.7579, 1.1401, 0.4945, 0.2122, 0.2613, 0.3984,
          14.0696,  10.6925,    0.0084, 8.5126, 2.1779, 1.1406, 0.4051, 0.3472, 0.5971,
        13.7417,    10.0652,    0.0084, 8.2977, 2.2172, 1.1204, 0.3929, 0.3583, 0.6143,
        11.7867,    15.1036,    0.0072, 8.3222, 1.9592, 0.9697, 0.3454, 0.3067, 0.5695),
      byrow=TRUE, ncol=9)
colnames(Y)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")

# mu.s values
mu.s<-matrix(c(6.97463333,  6.657444444,    10.560500000,   13.944108333,
          35.47273333,  13.221666667,   6.794500000,    4.149395833,
          0.00358000,   0.004533333,    0.008404167,    0.008891667,
          5.04523333,   7.128111111,    6.149708333,    6.932004167,
          1.19620000,   1.121222222,    2.086083333,    2.446995833,
          0.61910000,   0.485333333,    1.010833333,    0.880383333,
          0.19653333,   0.199222222,    0.491125000,    0.406904167,
          0.07113333,   0.312111111,    0.339916667,    0.298254167,
          0.45840000,   0.307222222,    0.616708333,    0.658416667),
         byrow=TRUE, ncol=4)
rownames(mu.s)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")
colnames(mu.s)<-c("CB","FD","RV","TS")

# Covariance matrix accounting for measurement error of control samples 
Sigma.inv<-structure(c(7.516889e-03,  4.719111e-04, 5.117778e-06,   4.896822e-03,   2.683444e-03,   5.368000e-04,   4.266667e-05,   -8.364444e-05,  4.630667e-04,
                  4.719111e-04, 8.097005e-05,   1.132174e-06,   8.006300e-04,   3.065903e-04,   4.421449e-05,   -4.869565e-06,  -5.677295e-06,  7.178551e-05,
                  5.117778e-06, 1.132174e-06,   9.927053e-08,   1.090671e-05,   4.145217e-06,   2.788406e-07,   -4.459903e-07,  -2.015942e-07,  9.478261e-07,
                  4.896822e-03, 8.006300e-04,   1.090671e-05,   9.350603e-03,   3.260499e-03,   4.947188e-04,   -3.441932e-05,  -4.135604e-05,  7.779478e-04,
                  2.683444e-03, 3.065903e-04,   4.145217e-06,   3.260499e-03,   1.390332e-03,   2.255014e-04,   -9.909179e-06,  -3.122995e-05,  2.906319e-04,
                  5.368000e-04, 4.421449e-05,   2.788406e-07,   4.947188e-04,   2.255014e-04,   7.714783e-05,   1.326377e-05,   -6.255072e-06,  4.525217e-05,
                  4.266667e-05, -4.869565e-06,  -4.459903e-07,  -3.441932e-05,  -9.909179e-06,  1.326377e-05,   1.484058e-05,   2.859903e-06,   -3.663768e-06,
                  -8.364444e-05,    -5.677295e-06,  -2.015942e-07,  -4.135604e-05,  -3.122995e-05,  -6.255072e-06,  2.859903e-06,   4.060386e-06,   -5.611594e-06,
                  4.630667e-04, 7.178551e-05,   9.478261e-07,   7.779478e-04,   2.906319e-04,   4.525217e-05,   -3.663768e-06,  -5.611594e-06,  7.074783e-05), dim=c(9,9))
Sigma.inv<-as.positive.definite(Sigma.inv)


# Sigma.s.inv values - a precision matrix
Sigma.s.inv<-structure(c(53.2149530065154, 7.78555807206446, -23455.5492262134, 
                     3.98226035437711, 181.147593576025, -554.144707627592, 810.959913329385, 
                     -222.049892052889, -676.160229178196, 7.78555807206446, 1.66555095380872, 
                     -2875.24335449319, 0.59057208979672, 35.2056764704436, -90.5513686454713, 
                     122.419754836375, -24.1464188840069, -110.985869289511, -23455.5492262134, 
                     -2875.24335449319, 26623099.2220065, -551.965187580738, -109728.586986558, 
                     276663.105767497, -386431.378216831, -82170.2904382779, 259038.414073531, 
                     3.98226035437711, 0.59057208979672, -551.965187580738, 0.918694682125983, 
                     9.31656679803296, -36.0723479510435, 38.501457515572, -39.2951967395694, 
                     -43.3516097843844, 181.147593576025, 35.2056764704436, -109728.586986558, 
                     9.31656679803296, 1046.08776695798, -2268.37491730957, 2936.52705816089, 
                     60.3905421599093, -2791.52606704426, -554.144707627592, -90.5513686454713, 
                     276663.105767497, -36.0723479510435, -2268.37491730957, 6336.11679874626, 
                     -8941.85134505843, 1263.96712724819, 7151.07272932156, 810.959913329385, 
                     122.419754836375, -386431.378216831, 38.501457515572, 2936.52705816089, 
                     -8941.85134505843, 14092.5591797243, -2587.63879760624, -10498.261488095, 
                     -222.049892052889, -24.1464188840069, -82170.2904382779, -39.2951967395694, 
                     60.3905421599093, 1263.96712724819, -2587.63879760624, 5072.16198995712, 
                     2677.23974494659, -676.160229178196, -110.985869289511, 259038.414073531, 
                     -43.3516097843844, -2791.52606704426, 7151.07272932156, -10498.261488095, 
                     2677.23974494659, 10377.5013426694, 39071331365.1695, -857206689.09299, 
                     13250725351142.8, 4293968882.52526, -138526471018.157, -346806050176.474, 
                     -185731096148.676, -40154997287.1756, -13832086879.5978, -857206689.09299, 
                     45035996.2737049, -291049796239.663, -94316253.2767046, 3042708991.03323, 
                     7617532441.71197, 4079550081.74548, 881997284.581011, 303819297.524648, 
                     13250725351142.8, -291049796239.663, 4499056063320832, 1457943322445.14, 
                     -47034282019304.7, -117752032879430, -63061801053651.4, -13633939080436.7, 
                     -4696447331898.41, 4293968882.52526, -94316253.2767046, 1457943322445.14, 
                     513861374.526114, -15241712287.5009, -38158180147.8331, -20435516111.352, 
                     -4418151354.45332, -1521909040.23616, -138526471018.157, 3042708991.03323, 
                     -47034282019304.7, -15241712287.5009, 491753442721.598, 1231009860799.41, 
                     659264193054.762, 142532685331.566, 49097861284.5838, -346806050176.474, 
                     7617532441.71197, -117752032879430, -38158180147.8331, 1231009860799.41, 
                     3081922920585.39, 1650491845606.43, 356835753178.233, 122918278499.906, 
                     -185731096148.676, 4079550081.74548, -63061801053651.4, -20435516111.352, 
                     659264193054.762, 1650491845606.43, 883961734443.138, 191102478025.129, 
                     65828570726.718, -40154997287.1756, 881997284.581011, -13633939080436.7, 
                     -4418151354.45332, 142532685331.566, 356835753178.233, 191102478025.129, 
                     41361275161.6651, 14232113705.0418, -13832086879.5978, 303819297.524648, 
                     -4696447331898.41, -1521909040.23616, 49097861284.5838, 122918278499.906, 
                     65828570726.718, 14232113705.0418, 4947128692.00388, 123.244158717727, 
                     15.1472733328015, -29469.8213424378, -162.933686694468, -830.16396214664, 
                     94.7853096969282, -68.8560246643179, 741.612274164014, 2264.45755928194, 
                     15.1472733328015, 5.35479274130578, -399.822164505876, -26.2908447667995, 
                     -48.3456870303433, -20.9894871349876, -7.0699744151633, 67.0416328725067, 
                     116.600590769204, -29469.8213424378, -399.822164505876, 25454078.4009846, 
                     34253.206150171, 238183.903396863, -33771.2583884048, 38523.5333162132, 
                     -288045.203177399, -1016390.47602203, -162.933686694468, -26.2908447667995, 
                     34253.206150171, 237.878595765458, 986.561801969833, -68.5166440497585, 
                     209.1509212556, -915.217874817308, -2932.38505457965, -830.16396214664, 
                     -48.3456870303433, 238183.903396863, 986.561801969833, 6669.29110276054, 
                     -1078.91990285319, 352.917095152379, -5609.73421176863, -17946.4520240348, 
                     94.7853096969282, -20.9894871349876, -33771.2583884048, -68.5166440497585, 
                     -1078.91990285319, 1052.57768971324, -913.546632402807, 484.410845136349, 
                     2431.79560931122, -68.8560246643179, -7.0699744151633, 38523.5333162132, 
                     209.1509212556, 352.917095152379, -913.546632402807, 2990.61062774375, 
                     525.038417043516, -3685.88997836577, 741.612274164014, 67.0416328725067, 
                     -288045.203177399, -915.217874817308, -5609.73421176863, 484.410845136349, 
                     525.038417043516, 7256.45767695138, 17207.7716622753, 2264.45755928194, 
                     116.600590769204, -1016390.47602203, -2932.38505457965, -17946.4520240348, 
                     2431.79560931122, -3685.88997836577, 17207.7716622753, 62661.2158837379, 
                     7.0546793436162, -0.0754213178285816, -2357.02875362113, 2.50749046482452, 
                     12.5166502527345, -86.6941792229212, 1.66923366612312, 7.93904545083449, 
                     -126.082128066627, -0.0754213178285816, 0.618360785191319, 1538.6589391227, 
                     -0.31799222497181, -1.26341570245756, 2.43556703274056, 18.6459686649981, 
                     -1.52020047556644, -1.33047266147014, -2357.02875362113, 1538.6589391227, 
                     7450231.40845542, -1467.47518565875, -8484.65408805661, 38482.6779430954, 
                     36677.1639056923, 5037.57654299884, -16285.3259202736, 2.50749046482452, 
                     -0.31799222497181, -1467.47518565875, 1.63035346006115, 1.13929112479727, 
                     -17.7642468677628, 2.98912100156106, 7.74881374361886, -54.902708360705, 
                     12.5166502527345, -1.26341570245756, -8484.65408805661, 1.13929112479727, 
                     106.706180072993, -244.571429668905, 24.7937851688084, -114.000942411739, 
                     -690.553652091766, -86.6941792229212, 2.43556703274056, 38482.6779430954, 
                     -17.7642468677628, -244.571429668905, 1534.27886445762, 549.368288309068, 
                     64.0901199201054, 980.291164452268, 1.66923366612312, 18.6459686649981, 
                     36677.1639056923, 2.98912100156106, 24.7937851688084, 549.368288309068, 
                     1781.88219182038, -107.322815612669, -1733.09607995999, 7.93904545083449, 
                     -1.52020047556644, 5037.57654299884, 7.74881374361886, -114.000942411739, 
                     64.0901199201054, -107.322815612669, 507.111649502059, 661.404359682356, 
                     -126.082128066627, -1.33047266147014, -16285.3259202736, -54.902708360705, 
                     -690.553652091766, 980.291164452268, -1733.09607995999, 661.404359682356, 
                     9442.32984843676), .Dim = c(9L, 9L, 4L))

#N
N<-nrow(Y)
#K
K<-4
#J
J<-9

myjagsfile<-"JAGSmodel.txt" # the above jags model
tV = t(ilrBase(D=K))
mydata <-list(N=N, J=J, K=K, Y=Y, mu.s=mu.s, Sigma.s.inv=Sigma.s.inv, R.Sigma=Sigma.inv, k.Sigma=J, tV=tV)
mypars <-c("Sigma.inv","p","kappa","phi", "s")
myiterations<-25000
myburnin<-12500
mythin<-175
set.seed(5000)
run.mcmc <- jags(data=mydata, inits=NULL, parameters.to.save=mypars, model.file=myjagsfile, n.chains=2, 
             n.iter=myiterations, n.burnin=myburnin, n.thin=mythin, DIC=TRUE)

Best Answer

The author of rjags has mentioned, that the Wishhart distribution can only be used, when it appears as a conjugate distribution in a strong sense:

http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/eab372de

There is a way around this though (I have to give credit to the authors of the R package bamdit though, thats where I first spotted the following trick:

model
{
for( i in 1 : n ) {
    m[i,1:2] ~ dmnorm(mu.0[1:2], sigma.inv[1:2, 1:2]) ## m is not the data
    w[i] ~ dgamma(nu.2, nu.2) T(0.1, 3)
    y[i, 1] <- mu[1] + m[i, 1] / sqrt(w[i]) ## y is the data
    y[i, 2] <- mu[2] + m[i, 2] / sqrt(w[i]) ## the scale enters here 
}
# Priors ...
mu[1] ~ dnorm(m.0[1], pre.mu[1])
mu[2] ~ dnorm(m.0[2], pre.mu[2])
mu.0[1] <- 0
mu.0[2] <- 0
# Weights distribution
nu.2 <- nu / 2
nu ~ dexp(nu.0) # prior for df
sigma.inv[1:2,1:2] ~ dwish(R[1:2,1:2], k)
sigma[1:2, 1:2] <- inverse(sigma.inv[1:2, 1:2])
}
Related Question