7
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ Are all of the steps Gibbs steps? If so, why not just write your own loop and forget about Rjags. $\endgroup$
    – user25658
    Commented Aug 2, 2013 at 14:51
  • $\begingroup$ I don't think the problem is necessarily Rjags, more that I am uncertain how best to represent measurement error within the model structure. As a relative newcomer to Bayesian statistics I find Rjags quite intuitive to use. $\endgroup$
    – Richard
    Commented Aug 2, 2013 at 15:23

1 Answer 1

3
$\begingroup$

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])
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.